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Abstract 


Orbital  altitudes  congested  with  spacecraft  and  debris  combined  with  recent  collisions 
have  all  but  negated  the  Big  Sky  Theory.  As  the  sheer  number  of  orbital  objects  to  track 
grows  unbounded  so  does  interest  in  prediction  methods  that  are  rapid  and  minimally 
computational.  Claimed  as  the  “other  solvable  solution,”  the  recently  completed  solution 
to  orbital  motion  about  the  earth,  based  on  Vinti’s  method  and  including  the  major  effects 
of  the  equatorial  bulge,  opens  up  the  prospect  of  much  more  accurate  analytical  models  for 
space  situational  awareness.  A  preliminary  examination  of  this  solution  is  presented.  A 
numerical  state  transition  matrix  is  found  using  Lagrange  partial  derivatives  to  implement 
a  nonlinear  least  squares  fitting  routine.  Orbit  fits  using  only  the  solvable  solution  for  non¬ 
circular,  non-equatorial  trajectories  less  than  60  degrees  inclination  are  on  the  order  of  a 
few  hundred  meters  with  projected,  average  error  growth  of  less  than  a  kilometer  per  day 
which  is  similar  to  the  expected  performance  of  the  Air  Force’s  method.  Also,  a  classical 
perturbations  approach  to  incorporate  the  dissipative  effects  of  air  drag  using  Hamiltonian 
action  and  angle  formulation  is  developed.  Predicted  drag  effects  are  97.5%  correct  after 
one  day  and  87%  correct  after  five  days  when  compared  to  an  integrated  truth.  Results  are 
validated  by  performing  a  similar  method  on  the  two  body  problem. 
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ORBIT  DETERMINATION  USING  VINTI’S  SOLUTION 


I.  Introduction 


1.1  Motivation 

On  12  February  2009,  an  active  US  communications  satellite,  Iridium-33,  and  an 
inactive  Russian  satellite,  Cosmos-2251,  collided  at  an  altitude  of  approximately  790 
km.  This  marked  the  first  known  accidental  collision  between  spacecraft  payloads.  The 
resulting  debris  was  estimated  to  range  between  1,000  and  2,000  objects  greater  than  10 
cm  and  between  60,000  and  120,000  objects  greater  than  1  cm  [1].  Later  that  same  year, 
the  USSTRATCOM  Commander  at  the  time,  General  Kevin  Chilton,  claimed  that  the  Big 
Space  Theory  came  to  a  close  with  this  seminal  event.  In  other  words,  no  longer  can  the 
space  community  claim  the  probability  of  collision  with  other  orbiting  objects,  due  to  the 
sheer  size  of  space,  is  low  enough  to  be  ignored  [2]. 

Due  to  orbital  changes  experienced  by  the  scattering  of  particles  upon  impact,  this 
significant  debris  cloud  has  been  dispersed  near  globally  and  now  congests  this  highly 
populated  altitude  regime.1  According  to  the  US  Air  Force,  as  of  1  March  2013,  it 
continued  to  track  2,160  pieces  from  this  collision  alone  [3].  This  number  is  approximately 
1,530  as  of  26  May  2015  [4], 

In  the  days  prior  to  the  collision,  a  close  approach  of  these  two  spacecraft 
was  predicted  by  SOCRATES.2  An  avoidance  maneuver  was  not  performed  prior  to 
the  collision,3  and,  with  perfect  hindsight,  it  could  be  surmised  this  was  completely 

'in  the  week  prior  to  the  incident,  it  was  reported  that  approximately  1,050  objects  would  come  within  5 
km  of  any  one  of  66  satellites  in  the  Iridium  constellation 

2  SOCRATES  is  a  public  service  to  the  international  satellite  community  provided  by  the  Center  for  Space 
Standards  &  Innovation  (CSSI)  -  www.celestrak.com/SOCRATES 

3  Only  Iridium-33  had  the  capability  for  maneuver 
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preventable.  However,  this  was  not  the  only  near  miss  for  that  week  and  not  even  the 
most  probable  or  closest  predicted.  The  miss  distance  predicted  for  that  conjunction  ranged 
between  117  m  and  1.812  km  across  reports,  while  the  closest  projected  approach  between 
several  other  crafts  averaged  between  30  and  50  m  [5]. 

It  would  seem  plausible  to  assume  that  if  a  craft  made  it  on  the  list  of  probable 
conjunctions,  an  avoidance  maneuver  should  be  planned  and  executed.  However, 
significant  mission  trade-offs  can  occur  with  even  the  smallest  bum.  Some  missions  require 
coming  off-line  for  a  period  of  time  causing  unneeded  outages.  Inevitably,  propellant  mass 
is  lost  that  is  ultimately  intended  to  maintain  the  payloads’  life  expectancy.  Therefore,  a 
serious  cost-benefit  analysis  is  required  for  each  possibility  of  a  close  call.  Planning  for  a 
maneuver  can  be  time-intensive  as  well,  requiring  several  hours  or  even  days.  What  makes 
this  decision  practical  is  an  accurate  estimate  of  both  crafts’  position  and  velocity  at  the 
future  time  in  question,  therefore  giving  a  realistic  probability  of  collision. 

Unfortunately,  using  the  data  available  to  determine  that  this  conjunction  was 
significant  in  the  midst  of  the  dozens  of  others  at  the  time,  is  “simply  not  possible” 
according  to  Kelso  [6].  It  is  surmised  the  amount  of  covariance  or  uncertainty  in  the 
position  estimate  becomes  greater  than  the  miss  distance  in  these  notices.  Statistical 
significance  of  this  information  is  lost  due  to  the  lack  of  actionable  information. 

An  increase  in  accuracy  of  current  satellite  state  estimates  combined  with  a  lack  of 
error  growth  over  time  would  revolutionize  conjunction  analysis.  The  vast  number  of 
false  positive  warnings  would  drastically  drop  and  statistics  would  become  reasonable. 
According  to  General  Robert  Kehler,  the  Commander  of  Air  Force  Space  Command  at 
the  time  of  the  collision,  10,000  conjunction  warning  messages  were  generated  by  the 
Joint  Space  Operations  Center  (JSPOC)  in  2012  but  only  75  maneuvers  were  performed 
in  response  [7].  Conservatively,  assume  those  75  instances  would  actually  have  resulted 
in  a  collision  if  no  action  were  taken.  This  means  the  warnings  issued  indicate  less  than 
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a  0.75%  chance  of  an  unplanned  rendezvous.  Realistically,  this  statistic  is  probably  much 
lower  but  the  calculated  probability  does  not  accompany  the  warning  [8].  If  accuracy  is 
improved,  these  messages  to  vehicle  operators  could  become  meaningful  with  enough  lead 
time  to  formulate  an  efficient  course  of  action. 

For  example,  consider  a  small  maneuver  of  1  cm/sec.  If  this  was  performed  one  hour 
before  conjunction,  it  would  result  in  a  position  change  of  36  m,  which  is  potentially  still  in 
harm’s  way.  However,  if  the  same  maneuver  is  accomplished  the  day  prior,  the  craft  would 
be  8.6  km  away  from  danger  at  that  time. 

1.1.1  Department  of  Defense  Concerns 
1.1. 1.1  Legacy  Analytical  Propagators 

The  US  Department  of  Defense  (DoD)  maintains  two  sets  of  mathematical  models 
to  calculate  and  predict  position/velocity  vectors  of  objects  in  low  earth  orbit  (LEO), 
Simplified  General  Perturbations-4  (SGP4)  and  Positions  and  Partials  with  Respect  to 
Time-3  (PPT3).  They  both  utilize  analytical  theory  and  are  originally  derived  and  modified 
from  Kozai  and  Brouwer’s  solution  (see  Chapter  II).  Using  these  programs  to  update  state 
vectors  for  all  foreseeable  objects  in  orbit  constitutes  maintaining  an  extensive  catalog  of 
satellites  and  orbital  debris  for  space  situational  awareness  purposes. 

With  the  exception  of  updates  to  account  for  more  objects  or  modernize  interface 
issues  and  computer  processes  within,  there  have  been  no  updates  to  the  underlying  theory 
in  over  50  years  [9].  Specifically,  Vetter  claims  the  last  “evolutionary  growth”  of  theory 
pertinent  to  SGP4  and  PPT3  occurred  in  1962  and  1964,  respectively.  To  put  this  in  military 
perspective,  air  superiority  workhorses  of  that  time  included  the  F-102  Delta  Dagger  and 
the  B-47  Stratojet.  Since  then,  significant  updates  to  orbital  solution  theory  have  been 
specific  to  special  perturbation  techniques  to  coincide  with  and  exploit  the  dawn  of  the 
computer  age  [10]. 
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1.1. 1.2  Accuracy  of  SGP4 

Although  the  actual  accuracy  of  SGP4  is  unknown  and  can  change  significantly  based 
on  several  factors  including  orbital  altitude  and  shape,  it  is  generally  known  to  be  not  very 
accurate.  Its  benefit  and  utility  have  long  been  advertised  as  a  first  cut  solution  to  give  to 
integrators.  Levit  labels  predictions  generated  from  SGP4  as  “not  sufficiently  accurate  to 
warrant  maneuvering  to  avoid  collision”  by  mentioning  it  creates  too  many  possibilities  of 
collisions  with  each  having  a  low  probability  [11].  This  could  be  likened  to  the  statement 
concerning  the  conjunction  assessment  lacking  statistical  meaning  (see  §1.1). 

Levit  goes  on  to  quantify  the  error  growth  of  a  new  method  as  100  m  per  day  and  says 
this  would  be  a  “10  fold”  improvement  upon  SGP4.  One  can  infer,  on  average,  that  SGP4 
can  be  expected  to  have  one  kilometer  of  error  growth  per  day.  Greene  introduces  a  method 
to  integrate  GPS  readings  into  SGP4  propagation  estimates  and  quantifies  unreasonable 
growth  rates  as  1  km  per  day  [12].  Dong  and  others  evaluate  an  unclassified  version  of 
SGP4  at  different  orbital  altitudes  with  different  parameters.  Solution  errors  for  near-earth 
objects  are  reported  to  be  on  the  order  of  a  few  hundred  meters  with  the  vast  majority  being 
under  10  km.  The  error  growth  over  the  sampled  orbits  is  determined  to  be  less  than  40  km 
over  3  days  for  LEO  satellites  and  only  1  day  for  those  with  elliptical  orbits  [13].  Vetter 
agrees  by  saying  the  analytical  theory  of  this,  and  similarly  PPT3,  are  only  accurate  to  a 
1-10  km  level  [10]. 

To  be  clear,  SGP4  does  not  calculate  miss  distance  for  specific  conjunctions.  Rather,  it 
provides  a  screening  of  all  active  satellites  for  possible  close  approaches.  If  certain  criteria 
are  met,  then  relevant  information  is  handed  over  to  more  accurate  systems  for  further 
analysis.  However,  these  other  algorithms  utilize  numerical  integration  and  take  much 
longer.  So,  if  the  first  cut  was  more  accurate,  time  spent  analyzing  meaningful  conjunctions 
would  be  minimized  by  avoiding  those  that  have  zero  chance.  System  time  could  be  opened 
up  for  other,  more  relevant  tasks. 
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1.1. 1.3  Orbital  Debris 


The  world’s  leading  orbital  tracking  entity,  the  US  Space  Surveillance  Network  (SSN), 
tracks  over  23,000  objects  in  LEO  larger  than  10  cm,  or  about  the  size  of  a  softball.  It  is 
estimated  there  are  upwards  of  500,000  objects  smaller  than  this  that  are  not  being  tracked. 
Of  all  the  objects  tracked,  only  about  5%  of  these,  or  about  1,300,  are  active  satellites  [14]. 
The  rest  are  considered  debris,  or  junk.  Currently,  to  maintain  the  catalog  of  over  23,000 
objects,  the  SSN  produces  over  200,000  sensor  taskings,  which  return  around  320,000 
observations  to  be  processed  every  day.  The  expected  activation  of  the  Space  Fence  at 
the  Kwajalein  Atoll  will  provide  the  ability  to  track  objects  as  small  as  2  cm.  It  is  clear 
that  the  amount  of  data  available  to  the  tracking  network  will  instantaneously  multiply 
astronomically. 

In  late  2014,  the  Orion  crew  capsule  orbited  the  earth  twice  and  then  experienced 
re-entry  as  part  of  an  experimental  test  flight  (Orion  EFT-1).  With  a  flight  duration  of 
almost  4.5  hours,  the  orbital  altitude  ranged  from  200  km  to  5,808  km.  Post-flight  analysis 
examined  the  recovered  portion4  of  the  craft’s  body  for  damage  and  found  25  regions  of 
interest.  After  further  investigation,  six  of  these  were  determined  to  be  potentially  due  to 
collisions  with  micro-meteroid  and  orbital  debris  (MMOD)  [15].  Although  the  size  of  these 
impact  craters  were  on  the  order  of  millimeters,  this  gives  reason  for  alarm.  A  flight  with 
such  short  duration  managed  to  collide  with  several  objects  that  are  too  small  to  track. 

1.1. 1.4  Debris  Growth 

Behind  risks  from  launch  and  deployment,  the  European  Space  Agency  (ESA)  claims 
space  debris  collision  ranks  as  the  third  highest  danger  to  many  missions  [16].  In  2008, 
Liou  and  Johnson  advertised  some  near-earth  orbital  regimes  had  reached  a  critical  density 
where  future  collisions  between  any  combination  of  active  satellite,  rocket  body,  or  debris 
would  be  inevitable  [17].  This  would  lead  to  a  creation  of  more  objects  kicking  off  a  chain 
4The  forward  bay  cover  was  jettisoned  during  parachute  deployment  and  sank  before  recovery 
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reaction  and  ultimately  create  an  unusable  orbital  region  consisting  of  a  debris  cloud.  This 
was  termed  the  Kessler  effect  in  1975  [18]. 

It  could  be  argued  that  the  realization  of  Liou’s  projection  occurred  with  the  Iridium- 
33/COSMOS-2251  collision.  To  provide  further  evidence,  on  22  January  2013  a  tiny 
Russian  satellite  named  BLITS5  collided  with  a  piece  of  debris  from  the  infamous  Chinese 
AS  AT  test  of  2007  [19].  These  are  not  the  only  collisions  to  have  occurred  in  recent  years, 
only  the  most  public. 

1.1. 1.5  Risk  Avoidance 

There  is  no  shortage  of  ideas  for  avoiding  the  risk  inherent  with  orbital  debris  growth. 
At  present,  methods  for  actual  debris  removal  from  orbit  are  just  that,  ideas.  The  Defense 
Advanced  Research  Projects  Agency  (DARPA)  has  initiated  the  Catcher’s  Mitt  study  to 
investigate  lowering  the  number  of  useless  objects  in  orbit.  Objectives  include  modeling  the 
orbital  debris  problem  coupled  with  its  expected  growth  and  exploring  technical  methods 
of  debris  remediation.  [20] 

In  support  of  this  study,  DARPA  commissioned  the  RAND  National  Defense  Research 
Institute  to  pose  the  problem  in  a  global,  strategic  light.  The  resulting  report,  Confronting 
Space  Debris,  compares  the  current  problem  to  other  globally  impacting  risks  such  as 
acid  rain,  airline  security,  and  oil  spills  to  name  a  few.  Characteristics  of  each  issue 
are  identified  within  categories  such  as  stakeholders,  blameworthy/affected  parties,  and 
mitigation/remediation  strategies. 

1.1.2  Recent  Findings 

The  National  Research  Council  of  the  National  Academies  published  a  study  named 
Continuing  Kepler’s  Quest:  Assessing  Air  Force  Space  Command’s  Astrodynamics 
Standards  in  2012  at  the  request  of  Air  Force  Space  Command  (AFSPC).  It  examined 
and  reviewed  current  standards  in  which  the  AF  uses  algorithms  and  computer  systems 

5Measuring  only  6.7  inches  in  diameter  and  weighing  a  mere  16.2  lbs,  BLITS  is  classified  as  a 
nanosatellite 
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to  maintain  awareness  of  thousands  of  orbiting  objects.  Among  several  suggestions  and 
findings  were  those  including  analytical  methods  and  the  command’s  understanding  of 
modern  dynamical  systems.  [8] 

Although  special  perturbations  has  widespread  popularity  for  modern  innovation  and 
accuracy  in  specialized  applications,  the  council  wanted  to  re-emphasize  the  importance 
of  updating  and  modernizing  general  perturbations  methods.  This  will  allow  continued 
interfacing  ability  with  international  satellite  users  while  decreasing  the  computational  load 
demanded  from  numerical  integrators  by  delivering  a  more  accurate  initial  solution. 

On  the  subject  of  the  organization’s  knowledge  of  modern  dynamical  theory,  the 
studies’  comments  were  scathing.  Not  only  did  the  committee  find  a  “striking  omission” 
in  modern  theory  employed,  they  were  not  convinced  the  experts  were  even  aware  of  such 
advancements.  It  was  reported  that  much  insight  into  the  behavior  of  dynamical  systems, 
which  is  critical  to  carrying  out  AFSPC’s  mission,  can  be  gained  by  no  longer  depending 
solely  upon  theory  untouched  since  the  1960s. 

1.2  Approach 

Over  50  years  ago,  Dr.  John  Pascal  Vinti  created  an  elegant  orbital  solution  using 
analytical  techniques  considered  advanced  even  by  today’s  standards  [21].  Due  to  the 
sophisticated  nature  of  his  techniques  and  a  lack  of  self  promotion,  this  field  and  the  DoD 
has  all  but  forgotten  his  effort  [22:1].  Wiesel  demonstrates  that  Vinti’s  solution  can  be  put 
through  modem  numerical  techniques  to  make  it  very  relevant  today  with  the  potential  to 
rival  current  methods  [23].  This  is  the  starting  point  for  the  current  work. 

First  and  foremost,  Wiesel’s  method  requires  a  state  transition  matrix.  Not  only 
is  it  required  by  least  squares  fitting,  it  is  needed  to  perform  meaningful  analysis  and 
development.  Although  analytically  finding  this  matrix  is  preferred,  the  author  proceeds 
with  a  version  developed  using  numerical  differentiation  techniques. 
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Once  local  orbital  motion  described  through  Vinti’s  solution  is  available  with  the 
state  transition  matrix,  it  can  be  analyzed  for  accuracy  and  viability.  A  least  squares  orbit 
fitting  routine  is  developed  and  employed  with  simulated  data  generated  using  a  numerical 
integration  truth  model.  Various  comparisons  are  first  made  using  Vinti’s  model  for  the 
earth’s  geopotential  (to  be  discussed  in  §2.5).  Observations  affected  by  the  full  geopotential 
(to  order  and  degree  20)  are  then  generated  and  fit  to  examine  initial  performance  that 
should  resemble  reality.  For  perspective,  this  performance  is  compared  to  that  of  an 
unclassified  version  of  SGP4  . 

As  a  demonstration  of  the  capability  inherent  in  Vinti’s  solution,  an  example  of  adding 
a  major  perturbation  is  presented.  Accounting  for  air  drag  is  demonstrated  using  the  action- 
angle  formulation  and  classical  perturbation  theory;  this  nonlinear  effect  is  sampled  and 
then  transformed  into  multiple  Fourier  series  that  can  be  analytically  integrated.  Having 
this  expression  provides  a  function  of  time  for  rapid  calculation  within  orbit  fitting.  For 
validation,  the  same  technique  is  performed  with  action-angle  variables  found  using  the 
two  body  problem  (TBP)  method. 

1.3  Problem  Statement 

The  recently  completed  solution  to  orbital  motion  about  the  earth,  based  on  Vinti’s 
method  and  including  the  major  effects  of  the  equatorial  bulge,  opens  up  the  prospect 
of  much  more  accurate  analytical  models  for  space  situational  awareness.  In  particular, 
prediction  of  conjunctions  absolutely  requires  accuracies  of  a  few  meters,  at  worst,  to 
eliminate  the  excessive  false  positive  rate  of  Space  Command  conjunction  messages.  The 
new  solution  will  be  evaluated  against  simulated  satellite  data  to  assess  its  accuracy.  With 
potential  upgrades  to  the  theory,  this  required  accuracy  for  collision  avoidance  could  be 
attained.  In  addition,  the  solution’s  utility  for  incorporating  additional  perturbations  such 
as  air  drag  is  examined. 
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1.4  Results 


When  the  newly  revised  orbital  solution  is  used  as  the  analytical  propagator  for  an  orbit 
fitting  routine,  results  are  promising.  For  orbits  only  perturbed  by  the  full  geopotential, 
errors  are  generally  on  the  order  of  hundreds  of  meters.  When  compared  to  an  unclassified 
version  of  SGP4,  results  were  only  worse  than  this  for  a  few  orbital  regimes.  Otherwise, 
Vinti  was  superior  or,  at  the  least,  very  similar.  Although  this  method  does  not  currently 
outperform  numerical  integration,  it  could  rival  it  after  further  perturbations  are  accounted 
for  in  future  works. 

The  problematic  exceptions  to  the  performance  described  above  are  when  inclination 
increases  past  65°.  Errors  rise  to  over  10  km  for  low  eccentricity  and  more  than  60  km  for 
high  eccentricity.  When  inclination  increases  past  82.5°  a  solution  is  no  longer  calculable 
due  to  a  numerical  flaw  in  the  integrals  deep  in  the  revised  solution.  Valid  results  return 
when  inclination  reaches  97.5°  and  mirrors  the  behavior  below  polar.  The  author  believes 
this  difficulty  is  due  to  the  negative  argument  of  perigee  rotation  rate  that  occurs  between 
the  two  critical  inclinations  of  63.4°  and  116.6".  Fixing  this  flaw  is  not  within  the  scope  of 
the  current  research. 

Also,  a  solution  is  not  available  for  eccentricity  and  inclination  of  exactly  zero  due 
to  the  root  finding  method  unique  to  the  new  theory  which  requires  a  definite  perigee 
and  northernmost  point.  However,  errors  did  not  grow  as  these  values  approached  zero. 
Eccentricity  and  inclination  as  low  as  0.001  and  0.01°,  respectively,  returned  satisfactory 
results. 

When  air  drag  is  accounted  for  through  a  novel  general  perturbations  approach,  the 
utility  of  Wiesel’s  solution  is  demonstrated.  Angle  perturbations  are  calculated  and  found 
to  be  97.5%  accurate  after  one  day  and  87%  after  five  days.  These  numbers  are  validated 
through  a  parallel  effort  using  TBP  variables. 
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1.5  Overview 


The  document  is  organized  into  seven  chapters.  This  first  chapter  has  provided  insight 
and  motivation  behind  what  makes  this  research  relevant.  Chapter  II  discusses  historical 
context  behind  the  development  of  analytical  orbital  solutions  as  it  pertains  to  space 
situational  awareness  and  catalog  maintenance.  Vinti’s  solution  and  modifications  made 
by  Wiesel  are  also  presented.  Chapter  III  introduces  a  method  for  applying  the  new  Vinti 
solution  to  an  orbit  fitting  routine.  A  numerical  state  transition  matrix  is  developed  and 
initial  fitting  results  are  examined.  Chapter  IV  summarizes  results  from  performing  orbit 
fitting  on  an  orbit  in  a  full  geopotential  using  the  Vinti  solution  as  modified  by  Wiesel  which 
uses  an  approximated  potential.  Chapter  V  presents  a  novel  approach  for  incorporating  air 
drag  perturbations  using  action-angle  variables.  This  approach  demonstrates  the  utility  in 
applying  Wiesel’s  solution  in  a  general  perturbations  context.  Chapter  VI  uses  the  approach 
presented  in  Chapter  V  but  with  TBP  variables  as  a  comparison.  Chapter  VII  provides 
conclusions  and  future  work  recommendations.  An  appendix  with  extra  graphs  exploring 
a  wider  range  of  orbital  parameters  is  included.  The  author  assumes  a  prior  knowledge 
of  Hamiltonian  mechanics  as  it  applies  to  orbital  dynamics  and  at  least  an  entry  level 
understanding  of  estimation  theory.  All  computer  code  is  developed  in  C++  with  double 
precision  accuracy  wherever  possible  in  units  of  nondimensional,  earth-centered  distance 
and  time  units. 
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II.  Background 


This  chapter  lays  the  groundwork  necessary  to  appreciate  using  the  Vinti  method  as 
a  modern  analytical  solution  for  earth-orbiting  objects,  be  it  spacecraft  or  debris.  A  brief 
history  of  orbit  determination,  as  it  applies  to  space  surveillance  and  tracking,  is  provided. 
Also,  a  few  classical  solution  methods  critical  to  that  history  are  discussed. 

2.1  Watching  the  Skies 

With  the  launch  of  Sputnik  in  1957,  the  military  mission  of  space  surveillance  was 
born.  Although  the  small  innocuous  object  only  transmitted  beeps  continuously,  it  implied 
that  an  adversary  on  the  other  side  of  the  globe  maintained  the  technological  capability  to 
reach  American  soil.  Therefore,  an  awareness  of  orbital  objects’  flight  paths  was  desired. 

The  country  needed  a  way  to  not  only  determine  the  location  of  an  enemy  spacecraft 
but  predict  where  it  would  be  in  the  foreseeable  future.  This  awareness  was  desired  by  the 
Navy  to  warn  battle  fleets  against  possible  reconnaissance.  Also,  if  a  satellite  could  pass 
overhead,  capabilities  definitely  existed  to  deliver  a  weapon  across  the  globe.  So,  as  the  Air 
Force  needed  to  provide  early  warning  for  such  incoming  threats,  object  differentiation  was 
critical.  It  was  paramount  that  they  had  the  capability  to  unequivocally  determine  between 
a  known  satellite  trajectory  and  an  inbound  warhead. 

To  this  end,  a  formal  effort  to  catalog  all  known  orbiting  objects  was  initiated  in  1959 
by  the  National  Space  Surveillance  Control  Center  (NSSCC)  at  Hanscom  Field  in  Bedford, 
MA.  Receiving  observations  by  way  of  teletype,  telephone,  mail,  and  personal  messenger, 
orbital  predictions  were  produced  using  relatively  primitive  methods,  then  distributed. 
These  predictions  included  estimated  times  of  passage  through  the  ascending  node  with 
corresponding  longitude  for  the  next  few  days.  Observation  sites  also  received  information 
regarding  where  to  look  in  the  sky  when  objects  flew  by  [24]. 
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Over  the  next  decade,  the  number  of  objects  in  orbit  rose  as  did  the  tide  of 
international  tensions.  Therefore,  the  community  desired  increased  accuracy  of  orbital 
predictions  through  more  sophisticated  methods.  Numerical  solutions  resulting  from 
special  perturbations  (to  be  discussed  in  §2.3.2)  provided  higher  accuracy  but  at  the  cost  of 
saturating  computer  systems  of  that  era. 

2.2  Evolution  of  Analytical  Theory  -  A  Snapshot 

Project  SPACETRACK  was  commissioned  to  seek  new  and  relevant  analytical 
solutions  describing  orbital  motion  to  support  catalog  maintenance.  In  the  same  edition  of 
Astronomical  Journal  in  November  of  1959,  Dirk  Brouwer  and  Yoshide  Kozai  published 
two  different  orbital  solutions  [25;  26].  Most  analytical  models  in  use  today  still  have  one 
of  these  two  original  formulations  as  their  foundation  [24]. 

Published  in  1959,  Brouwer’s  solution  to  orbital  motion  has  been  widely  used 
and  regarded  as  a  benchmark  of  artificial  satellite  theory.  He  performs  canonical 
transformations  to  simplify  the  equations  of  motion  using  the  canonical  Delauney  variables. 
In  terms  of  the  classical  orbital  elements,  these  are 


L  =  yfjla 

l  =  M 

G  =  LV  l-e2 

g  =  O) 

(2.1) 

H  -  G  cos  i 

h  =  Q. 

He  accomplishes  this  change  of  variables  using  a  geopotential  derived  from  a  spherical 
reference  frame,  shown  in  Figure  2.1.  Spherical  coordinates  are  related  to  the  cartesian 
frame  through 

x-r  cos  9  sin  0 

y  =  rsin0sin0  (2.2) 

2  =  rcos  (f) 
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Figure  2.1:  Spherical  Reference  Frame 


He  then  expands  the  disturbing  function  using  von  Ziepel’s  method.  This  provides 
short  and  long  period  terms  to  order  J2  and  secular  motion  to  order  The  Hamiltonian 
for  this  model  already  lacks  the  right  ascension  of  the  ascending  node,  Q,  but  his 
transformations  also  causes  argument  of  perigee,  co,  and  the  mean  anomaly,  M ,  to  drop  out. 
This  causes  the  conjugate  momenta,  L,  G,  H  to  be  constant.  Through  truncated  expansions, 
this  method  accounts  for  the  earth’s  asphericity,  J2  through  J5,  and  is  expressed  in  terms  of 
a  set  of  mean  elements  [27:688]. 

Kozai  generated  a  first  order  solution  that  has  rivaled  mainstay  orbital  solutions  to 
this  day.  He  uses  an  ad-hoc  averaging  technique  of  Lagrange’s  Variation  of  Parameters, 
which  is  inherently  fast  as  a  computing  method  [27:688].  In  1962,  he  went  on  to  derive  an 
extension  of  Brouwer’s  original  solution  [28]. 

Also  in  1959,  Garfinkel  attempted  a  similar  method  to  Brouwer  [29].  Using  von 
Ziepel’s  method  with  Delauney  variables,  the  most  noticeable  departure  in  theory  is  the 
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form  of  the  geopotential  used.  However,  Garfinkel  claimed  his  results  were  similar  to  those 
of  Brouwer. 

Over  the  next  decade,  development  of  analytical  orbital  theory  continued  mostly  in 
the  context  of  modifying  and  adding  to  the  original  solutions  of  Brouwer  and  Kozai.  It  was 
in  these  solutions  that  the  US  Air  Force’s  analytical  orbital  solution  software,  SGP,  was 
rooted  with  its  genesis  in  the  early  1960s  [24]. 

Lyddane  significantly  improved  Brouwer’s  solution  by  implementing  Poincare 
variables  with  which  to  solve  the  system  [30].  Before,  singularities  existed  at  low 
inclination  and  eccentricity  and  at  the  critical  inclination.  This  change  of  variables 
preserved  the  benefits  of  the  solution  while  avoiding  problematic  results.  Similar  to 
Poincare’s  approach,  Deprit  addressed  problems  at  the  critical  inclination  [31]. 

Lyddane’s  modification  of  Brouwer’s  solution  became  the  defining  characteristic 
that  set  apart  the  US  Navy’s  analytical  propagator,  PPT,  or  Positions  and  Partials  with 
respect  to  Time.  The  US  Air  Force’s  SGP  continued  with  a  modified  blend  of  Kozai 
and  Brouwer’s  theory  until  later  adopting  Lyddane’s  modification  in  SGP4  [24].  Brouwer 
and  Hori  attempted  to  incorporate  drag  effects  by  using  a  static  exponential  representation 
for  atmospheric  density  with  a  constant  scale  height  [32].  This  complex  model  was  so 
extensive  for  computers  at  the  time  that  it  was  operationally  infeasible.  Simplified  treatment 
of  the  density  profile  allowed  improvements  first  by  Lane,  then  by  Lane  and  Cranford 
[33;  34].  These  improvements  as  well  as  Lyddane  and  Deprit’s  modification  were  key  to 
the  development  and  implementation  of  SGP4. 

Significant  improvement  in  orbital  theory  and  operations  occurred  in  this  turbulent 
decade  of  the  1960s.  However,  in  these  early  years  of  space  tracking,  solutions  were 
accurate  enough  if  they  simply  provided  acceptable  look  angles  and  flyover  times. 
Conjunctions  or  close  flybys  of  other  spacecraft  or  objects  were  almost  inconceivable  due  to 
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the  Big  Space  Theory,  a  derivative  of  the  Big  Sky  Theory.  Concerning  these  improvements 

specific  to  Brouwer’s  theory,  Vallado  says: 

All  these  developments  give  Brouwer’s  method  the  appearance  of  being  a 
superb  analytical  theory,  but  it  has  several  significant  shortcomings.  It’s  fine 
for  general  use  in  applications  needing  limited  accuracy.  But  as  computational 
power  increases  and  satellite  systems  demand  more  accuracy,  its  effectiveness 
diminishes  -  it’s  still  an  analytical  series  approximation  accounting  for  a 
large  but  incomplete  subset  of  dominant  perturbations,  so  it’s  only  moderately 
accurate.  The  U.S.  Air  Force  and  Navy  operate  similar  versions  of  this  theory. 
[27:690] 

So,  in  this  modem  context,  the  community  should  not  remain  complacent  with  the 
confidence  in  a  traditional  method  merely  because  it  has  been  used  for  decades. 

In  the  past,  if  more  accuracy  was  needed  for  a  certain  orbit,  analytical  theory  was 
put  aside  after  a  first  cut  solution  was  obtained.  The  numbers  were  then  handed  over 
to  numerical  integrators.  This  is  still  true  today,  however,  in  the  past  the  solution  was 
more  accurate  at  the  expense  of  computing  time  for  the  rest  of  the  afternoon.  Although 
computing  time  necessary  for  this  alternate  approach  has  reduced  significantly,  rapid 
analytical  solutions  still  have  value.  This  is  especially  true  for  those  that  are  actually 
accurate  enough  to  be  meaningful.  Using  these  two  different  methods  for  orbital  prediction 
will  be  discussed  in  the  next  section. 

2.3  Orbit  Determination  Methods 

Determining  the  unique  elements  of  a  satellite’s  orbit  is  only  as  useful  as  how  well 
those  elements  can  be  used  to  predict  where  that  craft  will  be  at  a  specific  time  in  the 
future.  The  accuracy  of  such  a  forecast  depends  upon  the  model  used  in  such  a  calculation. 
Predominantly,  there  are  three  types  of  orbit  fitting/determination:  general  perturbations, 
special  perturbations,  and  semianalytical  techniques.  The  latter  combines  aspects  of  the 
first  two  and  will  not  be  reviewed  in  this  work  due  to  the  wide  variety  of  techniques  used 
and  lack  of  extensive  documentation.  Only  the  fundamental  premises  behind  the  first  two 
approaches  will  be  presented. 
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2.3.1  General  Perturbations 


Also  known  as  analytical  theory,  Vallado  explains  that  “General  Perturbations 
techniques  replace  the  original  equations  of  motion  with  an  analytical  approximation  that 
captures  the  essential  character  of  the  motion  over  some  limited  time  interval  and  which 
also  permits  analytical  integration”  [27:609].  So,  instead  of  describing  the  unique  motion 
with  exact  solutions,  successive  approximations  are  made  using  various  techniques.  These 
expressions  are  generally  valid  for  all  initial  conditions  but  only  for  a  limited  amount 
of  time.  Not  unique  to  orbital  mechanics,  many  disciplines  utilize  this  methodology  in 
obtaining  useful  physical  descriptions  of  their  systems  of  choice.  The  solutions  of  Kozai, 
Brouwer,  and  Garfinkel  fall  into  this  category. 

However,  expanding  the  expressions  analytically  can  quickly  get  overly  complex. 
More  often  than  not,  only  second  order  solutions  are  provided  due  to  these  difficulties. 
Neglecting  higher-order  effects  obviously  degrades  the  accuracy  but  it  makes  a  solution 
attainable.  Also,  describing  the  motion  in  such  a  way  can  provide  valuable  insight  into 
the  inner  workings  of  the  dynamics  to  aid  in  stability  or  ancillary  studies.  Although  the 
derivations  are  arduous,  computing  time  for  a  solution  is  minimal.  This  analytical  method 
was  commonplace  before  modern  computing  techniques  opened  the  possibility  of  rapid 
numerical  integration,  bringing  special  perturbations  to  the  forefront. 

2.3.2  Special  Perturbations 

In  the  modem  age  of  computing,  special  perturbations  has  been  the  leader  of  accurate 
OD  schemes  and  orbital  propagation.  This  involves  capturing  perturbations  as  accelerations 
in  accurate  equations  of  motion  that  can  be  numerically  integrated  using  a  wide  range  of 
methods  and  schemes.  A  specific  orbit  is  defined  by  a  set  of  initial  conditions  and  then  is 
integrated  until  a  time  of  interest  is  reached.  But,  as  the  name  would  indicate,  a  particular 
solution  is  unique,  or  special,  and  cannot  be  used  for  analyzing  other  orbits;  the  process 
must  be  re-initiated  for  a  new  set  of  initial  conditions.  This  is  unlike  general  perturbations  , 
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where  an  analytical  solution  is  valid  for  any  set  of  initial  conditions.  Also,  due  to  truncation 
and  roundoff  errors  inherent  in  the  integration,  errors  build  up  as  the  square  root  of  the 
number  of  calculations  made  [35:2]. 

Coffey  has  introduced  a  way  for  the  Navy  to  perform  catalog  maintenance  using 
special  perturbations  [36].  He  contends  that  accuracy  will  increase  while  mission 
requirements  will  be  maintained.  He  uses  a  parallelized  computing  approach  and  says 
to  “just  add  more  [computer]  processors”  when  the  immense  volume  of  tracked  objects 
saturates  the  system.  This  statement  about  tracking  all  the  objects  with  SP  was  made  when 
the  cataloged  items  was  on  the  order  of  20, 000.  If  the  nation  begins  tracking  objects  as 
small  as  2  cm  after  bringing  a  new  Space  Fence  online,  the  catalog  would  suddenly  grow 
to  over  half  a  million  and  instantly  multiply  processing  time  for  any  system  attempting  to 
track  them  all  [14]. 

2.4  Space  Surveillance  and  Tracking 

Analytical  propagators  help  maintain  situational  awareness  over  all  objects  in  the 
catalog.  They  also  provide  various  tools  and  information  to  the  satellite  community, 
including  look  angles6  and  initial  conjunction  analyses.  When  increased  accuracy  for  a 
few  objects  is  needed,  the  trajectories  are  analyzed  further  by  integrators. 

In  keeping  a  database  of  thousands  of  objects,  general  perturbations  wins  out  in 
keeping  up  with  them  all  because  only  one  computer  call  is  necessary  to  provide  state 
information  for  a  given  time  in  the  future.  Special  perturbations  would  have  to  perform 
calculations  at  each  time  step  from  epoch  until  desired  time.  Although  the  estimate  would 
probably  be  more  accurate,  even  today’s  computers  could  potentially  be  overburdened 
performing  such  operations  on  thousands  of  satellites  simultaneously. 

The  reduced  complexity  in  computing  combined  with  the  ability  to  quickly  approx¬ 
imate  a  state  vector  sometime  in  the  future  surely  makes  analytical  theory  attractive  for 
6This  tells  the  user  where  to  point  the  antenna  for  uplink/downlink  purposes 
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space  situational  awareness  as  it  applies  to  catalog  maintenance.  The  current  research  ex¬ 
amines  a  new  solution  that  maintains  computing  benefits  of  general  perturbations  yet  with 
potential  to  approach  numerical  integration  accuracy. 

2.5  Vinti’s  Solution 

Return  to  that  pivotal  point  in  history  when  the  US  was  racing  to  modify  current 
methods  to  maintain  situational  awareness  in  the  race  for  space.  The  preponderance  of 
analytical  solutions  utilized  spherical  coordinates  as  discussed  earlier.  Using  this  reference 
frame,  the  Hamiltonian  for  two  body  motion  can  be  separated  and  give  a  starting  point 
from  which  to  begin  adding  perturbations.  However,  this  approach  only  results  in  one 
fundamental  frequency  from  which  to  begin,  termed  fully  degenerate. 

Going  almost  unnoticed  in  the  shadow  of  the  likes  of  Brouwer  and  Kozai,  John  Pascal 
Vinti  posed  a  unique  departure  in  theory.  Leveraging  his  physics  background,  he  introduced 
an  elegant  oblate-spheroidal  coordinate  system  to  describe  the  orbital  motion  as  perturbed 
by  the  earth’s  zonal  effects  [21].  Illustrated  in  Figure  2.2,  the  spheroidal  coordinates  are 
related  to  the  physical  ones  through  the  expressions 

x  =  V(P2  +  c2)0  -  H2) cos  0 

y  =  y]{p2  +  c2)(l  -  if)  sin  (p  (2.3) 

z  =  pn 

where  p  determines  the  size  of  the  spheroid  and  is  approximately  the  radius  and  c  is  a 
constant  that  determines  the  shape.  For  earth,  c  is  set  equal  to  J2  times  the  radius  of  the 
earth  squared.  The  right  ascension  of  the  satellite  is  (p  and  //  determines  the  shape  of  the 
confocal  hyperboloids  of  revolution. 

Due  to  the  orthogonality  introduced  between  the  reference  frame  and  the  effects  due 
to  zonal  influence,  developing  the  gravitational  potential  using  this  reference  frame  exactly 
accounts  for  J2  and  almost  75%  of  J 4.  Then,  performing  a  change  of  variables  within  this 
frame  leads  to  separability  of  the  perturbed  Hamiltonian- Jacobi  equation.  Separating  this 
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equation  was  Brouwer  and  Kozai’s  motivation;  however,  they  accomplished  the  needed 
change  of  variables  only  through  expansions,  which  limit  accuracy.  The  reader  is  directed 
to  Appendix  A  for  the  Hamiltonian  expressions  resulting  from  Vinti’s  method. 


Although  the  solution  was  unique  and  elegant,  difficulty  arose  with  the  inherent 
formulation  and  use  of  elliptical  integrals  in  manipulating  such  a  solution.  These  are  found 
in  the  resulting  generating  function 


W  =  a^A : 


y/F(p,  a) 

r.2  i  r>2  dP 


y/G(T],  a) 

— ; - —dri 

1  —  77- 


(2.4) 


cz  +  pz 

where  the  integration  contours  p~  and  rj~  indicate  integrating  once  around  the  orbit 
beginning  at  perigee  and  the  southernmost  point,  respectively.  Also,  F  and  G  are  quartics 
of  the  form 

F(p,  a)  =  c2al  +  (c2  +  p2)(-a\  +  2 pp  +  2 (aq  +  co^a^)p2)  (2.5) 

G(?j,  a)  =  - a 2  +  (1  -  rf)(a\  +  2(ai  +  m®cr3)cV)  (2-6) 
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with  at  being  the  separation  momenta.  At  the  time,  computing  methods  were  inadequate 
at  best  for  handling  such  expressions.  Solving  these  elliptical  integrals  analytically 
encouraged  the  community  to  adopt  the  other  seemingly  viable  and  more  usable  solutions 
present  at  the  time.  [23] 

This  lack  of  interest  did  not  discourage  Vinti  from  continuing  his  attempts  to  solve 
and  improve  upon  his  solution.  In  1961,  he  published  the  derivation  of  the  theory  and 
then  somewhat  of  a  recipe  book  for  computing  non-equatorial  intermediary  orbits  using  his 
solution  [37;  38].  In  1962,  he  furnished  a  corollary  to  the  previous  papers  by  addressing 
equatorial  orbits  [39].  In  1963  he  introduced  a  method  to  account  for  zonal  perturbations 
to  his  solution  [40].  These  modifications  were  superseded  in  1966  when  he  adjusted  his 
original  solution  to  account  for  the  third  zonal  harmonic  [41].  A  subsequent  improvement 
was  released  in  1969  but  was  structured  very  similarly  to  the  recipe  book  of  1961  [42]. 
Using  the  new  solution,  he  goes  on  to  demonstrate  how  drag  can  be  implemented  using  his 
unique  variables  [43]. 

The  commonality  in  all  the  works  presented  here  lies  in  the  oblate  spheroidal  solution 
method.  Accomplishing  the  change  of  variables  required  to  separate  the  Hamilton- Jacobi 
equation  is  accomplished  directly  without  expansions.  Not  only  is  this  unique  for  orbital 
mechanics,  it  is  remarkable  considering  that  the  full  effect  of  J2  and  the  majority  of 
J4  are  included  in  the  potential.  However,  the  resulting  equations  of  motion  involving 
the  infamous  elliptical  integrals  are  solved  by  expansions  in  J2,  thus  reducing  potential 
accuracy. 

Recently  edited  in  1998,  Vinti’s  published  textbook  furnishes  the  latest  in  attempts  at 
a  closed  form  solution  to  his  problem  posed  almost  40  years  prior  [22:99-103].  However, 
these  are  still  accomplished  using  expansions  and  thus  introduce  inaccuracies  and  unneeded 
complexities  as  opposed  to  solving  the  problem  directly;  doing  so  is  demonstrated  by 
Wiesel  (see  §2.7). 
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2.6  Works  Acknowledging  and  Improving  Vinti’s  Method 

Since  its  publishing  in  1959,  the  majority  of  authors  who  cite  Vinti’s  solution 
acknowledge  the  elegance  in  such  a  formulation  but  quickly  take  a  different  tack  due  to 
the  difficulties  inherent  in  actually  solving  it.  For  instance,  Barrar,  in  1961,  claims  that 
“...Vinti’s  solution  is  undoubtedly  an  excellent  one.  However,  Vinti  only  reduced  the 
solution  to  quadratures.  His  result  is  four  elliptic  integrals  of  the  third  kind.”  He  then 
proceeds  to  present  a  different  method  that  ends  up  resembling  Garfmkel’s  solution  [44]. 
More  recently,  in  2012,  Morrison  even  goes  as  far  as  to  say  that  these  resulting  expressions 
result  in  “horrendous  analysis,  which  we  have  pledged  to  avoid.”  [45:232] 

A  much  smaller  population  are  those  who  acknowledge  the  potential  of  using  such  a 
solution  in  an  operational  application.  For  example,  Izsak  starts  with  Vinti’s  “remarkable 
approximation”  to  earth’s  geopotential  and  attempts  to  perform  perturbation  theory  to 
account  for  the  difference  between  this  and  reality.  The  magnitude  of  this  difference  is 
on  the  order  of  His  solution  claims  to  account  for  J2  through  J$,  much  like  Brouwer’s 
[46], 

In  1964,  Allen  and  Knolle  take  Izsak’s  formulation  and  develop  a  method  of  orbit 
fitting  by  calculating  differential  correction  coefficients  for  least  squares  [47].  This 
approach  is  very  similar  to  what  will  be  presented  in  Chapter  III.  In  1967,  Walden  and 
Watson  also  applies  differential  corrections  to  this  problem  but  focuses  on  Vinti’s  later 
modification  in  1966  which  includes  the  third  zonal  harmonic  [48]. 

First  in  1962  for  the  original  solution,  then  again  in  1966  for  the  inclusion  of  the  third 
harmonic,  Bonavito  provides  computational  procedures  for  implementing  Vinti’s  theory  for 
computing  an  accurate  reference  or  intermediary  orbit  [49;  50].  Then  joining  with  others 
in  1969,  he  compares  the  theories  of  Brouwer  and  Vinti  in  accuracy  and  speed.  To  be  clear, 
this  Brouwer  theory  is  the  original  one  from  1959  and  not  the  one  blended  with  Kozai  and 
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others  to  form  SGP4.  Bonavito’s  findings  show  that  Brouwer’s  theory  is  “considerably  less 
accurate”  than  that  of  Vinti  [51]. 

Similar  to  Bonavito’s  comparison,  Gordon  et  al  examine  differences  between  Brouwer 
and  Vinti  theory.  However,  they  also  consider  Brouwer’s  solution  as  modified  by  Lyddane. 
They  present  results  that  confirm  Vinti’s  superiority  except  for  cross-track  errors.  They 
explain  this  is  believed  to  result  from  the  J\  approximation  for  J4,  which  is  about  70%  of 
the  value.  The  work  then  presents  a  technique  to  perform  a  first-order  general  perturbations 
approach  to  account  for  this  difference.  After  this  change  is  made,  results  significantly 
improve.  Cross-track  errors  are  improved  by  more  than  96%  and  in-track  errors  were 
reduced  by  90%  in  certain  orbits  [52]. 

One  of  the  few  remaining  modern  day  proponents  of  the  Vinti  theory  is  Gim  Der.  In 
1996,  he  provided  derivations  and  comparisons  for  Keplerian,  Vinti,  and  numerical  state 
transition  matrices  [53;  54],  Joining  the  likes  of  Bonavito,  he  edits  Vinti’s  textbook  in  1998 
and  provides  an  introduction  applauding  the  “brilliant  effort”  of  Vinti  [22:  viii] .  Founding 
DerAstrodynamics  in  2002,  Der  “...provid[es]  innovative  and  specialized  Astrodynamics 
algorithms  for  Space  Situational  Awareness  applications”  [55].  The  theory  of  his 
algorithms  and  approach  is  rooted  in  Vinti’s  method. 

With  the  exception  of  Der  and  Wiesel,  Vinti’s  solution  has  more  or  less  been  forgotten. 
All  the  related  efforts  mentioned  above  use  Vinti’s  method  and  the  improvements  thereafter 
as  the  basis  for  their  research.  The  focus  of  the  current  research  is  to  examine  Wiesel’s 
approach  that  results  from  modifying  Vinti’s  original  theory  in  1959.  His  solution  is,  in 
every  respect,  a  new  one  and  should  be  treated  as  such. 

2.7  Vinti’s  Problem  As  Modified  By  Wiesel 

Wiesel  takes  Vinti’s  problem  that  has  laid  dormant  for  several  decades  and  examines 
it  through  the  lens  of  modem  computing.  Instead  of  posing  that  the  elliptic  integrals  be 
performed  analytically  or  through  expansions,  he  offers  they  can  be  solved  numerically. 
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Before  giving  the  task  to  the  computer,  he  manipulates  the  solution  forms  to  avoid 
singularities  and  to  provide  more  usable  structures  upon  which  to  perform  further  beneficial 
calculations  [23]. 

In  solving  Vinti’s  problem  directly,  Wiesel  avoids  expansions  in  J2  that  limited  Vinti’s 
potential  accuracy  when  he  attempted  to  solve  his  integrals  analytically.  Wiesel  performs 
a  change  of  variables  to  remove  square  root  singularities  in  the  elliptic  expressions  thus 
providing  an  avenue  to  furnish  an  obtainable  numerical  solution.  Considering  a  quartic  can 
be  factored  into  its  leading  coefficient  and  roots  as 

F(p,  a)  =  A(p  -  pi)(p  -  p2)(p+  -  p)(p  -  p~)  (2.7) 


and  introducing  the  change  of  variables  illustrated  in  Figure  2.3  as 


P  = 


P+  +P 
2 


cos  Ep 


(2.8) 


Figure  2.3:  Circumscribing  A  Circle  Between  Integration  Limits 


Wiesel  shows  the  square  root  factor  of  the  quartic  above  (Equation  2.5)  becomes 

y/F(p,a)  =  -\jApip  —  pi)(p  —  p2)  \p+  -p~ |  |sin£p|/2  (2.9) 
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where  p\  and  p2  are  the  roots  that  do  not  bracket  the  real  trajectory  and  Ap  is  the  coefficient 
of  p4.  He  repeats  this  process  for  G(rj,  a).  In  doing  so,  he  generates  two  eccentric  anomaly¬ 
like  expressions  for  Ep  and  Et]  requiring  simultaneous  computation.  This  new  root  finding 
method  creates  numerical  difficulties  at  exactly  equatorial  (lacking  a  definite  southernmost 
point)  and  circular  (lacking  a  definite  perigee  point)  orbits  as  well  as  near  polar  (argument 
of  perigee  rate  is  negative)  orbits.  Details  of  this  are  discussed  in  Chapter  IV.  These 
challenges  could  potentially  be  circumvented  in  the  future  through  somewhat  of  a  patch 
by  reverting  to  Vinti’s  expansions  in  the  neighborhood  of  such  orbits.  Providing  such  a  fix 
is  outside  the  scope  of  the  current  research  but  is  recommended  for  future  work. 

Wiesel’s  process  creates  a  new  solvable  problem  from  which  to  begin  perturbations. 
He  claims  that  this  is  “the  other  solvable  problem.”  In  comparison  to  the  widely  used  two 
body  problem,  perturbations  can  begin  on  the  order  of  10-5  instead  of  10-3. 

Not  only  is  Vinti’s  solution  separable  but  it  can  provide  action-angle  variables  as 
promised  by  Hamilton-Jacobi  theory.  Having  the  generating  function  in  hand,  Wiesel 
demonstrates  they  can  be  found  through  performing  a  contour  integral  as 


li 


i  X  ,  i  r  aw 

—  ( )  p4q:  =  —  <T>  - 

2^r  Jr, 


-dcfr 


(2.10) 


In  Jj-  dqt 

where  T,  is  a  contour  around  the  torus  encircling  only  one  coordinate.  To  visualize  this, 
consider  a  two  dimensional  system  represented  in  phase  space.  Figure  2.4  shows  this 
becomes  a  three  dimensional  torus  in  four  space  and  how  the  integral  goes  around  each 
coordinate  separately. 

For  the  quartic  examined  above,  this  expression  looks  like 


In 


1  f  7 F(p ,  a ) 


■2  nfp 


c2  +  p2 


dp 


(2.11) 


Separability  of  a  system  where  time  is  independent  indicates  that  each  integral  of  the 
motion  only  depend  on  one  variable  and  are  therefore,  uncoupled.  Having  access  to 
the  system’s  action-angle  formulation  provides  valuable  insight  into  periodic  systems  by 
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furnishing  basis  frequencies.  Obtaining  these  frequencies  using  Hamilton- Jacobi  theory 
has  been  termed  Delauney’s  method  [56:361-372].  Thus,  the  Delauney  variables  of  orbital 
motion  are  the  result  of  the  two  body  problem  being  solved  through  such  a  process. 


Figure  2.4:  Contour  Example:  2-torus  in  4-space 


Access  to  these  valuable  quantities  is  key  to  describing  orbits  on  KAM  tori.  A  brief 
introduction  to  KAM  theory  is  offered  in  §2.8.  Utilizing  this  theory  will  not  be  the  emphasis 
of  the  current  research,  but  using  action-angles  in  a  classical  perturbation  approach  will  be 
demonstrated  for  air  drag  in  Chapters  V  and  VI. 

Duffy  explores  orbital  resonances  within  Vinti’s  solution  [57].  Using  truth  orbits  for 
GPS  satellites  that  sit  on  the  2:1  resonance,  she  examines  stable,  librational,  and  chaotic 
behaviors  over  different  ten  year  trajectories.  To  avoid  small  divisors  within  the  Vinti 
solution  for  these  resonances,  she  presents  a  new  transformation  of  variables. 

A  brief  note  on  reference  frames  is  in  order.  Wiesel  performs  an  extra  transformation 
in  order  to  introduce  a  rotating  reference  frame.  This  sets  the  stage  for  using  Vinti’s  solution 
to  express  orbits  as  lying  on  KAM  tori.  Since  this  is  not  the  primary  focus  of  the  current 
research,  a>@  will  be  set  to  zero,  thus  providing  an  inertial  frame  [23]. 
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2.8  KAM  Theory 

For  perspective,  a  brief  explanation  and  history  of  KAM  theory  is  provided.  First 
illustrated  by  Kolmogorov  then  later  proved  by  Arnold  and  Moser  (thus  providing  the 
theory’s  namesake),  the  KAM  theory  has  mainly  laid  dormant  in  regards  to  orbital  problems 
[58-61].  Introduced  over  50  years  ago,  it  offers  that  solutions  for  most  lightly  perturbed 
Hamiltonian  systems  lie  geometrically  on  tori  in  phase  space. 

In  the  last  decade,  KAM  theory  has  enjoyed  a  revival  in  being  applied  to  modern 
orbital  dynamics  problems.  The  work  of  Celletti  and  Chierchia  applies  it  to  the  Planar, 
Circular,  Restricted  Three  Body  system  of  the  Sun,  Jupiter,  and  the  asteroid  Victoria  [62]. 
The  remaining  related  works  are  predominantly  by  Wiesel  and  his  graduate  and  post¬ 
graduate  students  at  AFIT.  These  examine  the  theory  as  it  applies  to  artificial  satellites 
and  comprise  the  majority  of  the  literature  with  this  purpose. 

By  taking  the  theory  further,  Wiesel  suggests  that  earth-orbiting  satellites  lie  on  KAM 
tori  [63].  This  is  a  new  and  promising  lens  through  which  to  view  artificial  satellite  theory. 
The  work  by  Little  takes  Wiesel’s  approach  and  attempts  to  compute  tori  on  which  actual 
satellites  lie  so  that  accurate  positions  are  known  [64].  Visher  provides  an  independent 
verification  of  using  tori  for  orbit  propagation  by  comparing  results  to  Systems  Tool  Kit 
(STK),  as  a  truth  model  [65].  Further  examination  of  various  methods  to  construct  a  torus 
through  new  spectral  methods  is  preformed  by  Wiesel  and  Bordner  [66;  67]. 

Craft  explores  the  utility  and  viability  in  planning  for  satellite  constellations  to  exist  on 
the  same  torus,  only  separated  by  angle  displacements.  The  benefit  here  is  the  near  constant 
formation  separation  reducing  fuel  allowance  for  constellation  station-keeping[68].  Yates 
explores  the  compensation  errors  in  reference  tori  for  effects  due  to  air  drag  and  third  body 
perturbations  [69].  Hagen  goes  on  to  examine  these  effects  as  it  applies  to  KAM  theory 
and  shows  that  the  theory  still  holds  in  the  presence  of  these  light  perturbations  [70].  Dunk 
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applies  lessons  and  methods  learned  earlier  to  the  unique  environment  of  highly  eccentric 
orbits  [71]. 

Frey  demonstrated  extracting  KAM  torus  basis  freqencies  from  SGP4  and  two- 
line  elements  sets  (TLE).  However,  this  method  is  problematic  at  low  eccentricity  [72]. 
Using  Wiesel’s  new  theory  for  nearly  circular  orbits  [73],  Abay  performs  orbit  fitting  by 
constructing  KAM  tori  based  on  TLEs,  of  orbits  with  low  eccentricity.  When  compared  to 
SGP4  fits  for  the  same  orbits,  results  are  up  to  5  times  better  [74]. 

This  review  of  literature  relating  to  and  applying  KAM  theory  to  artificial  satellite 
motion  is  not  all  inclusive.  It  does,  however,  highlight  the  emergence  of  modern 
applicability  and  practicality  to  a  theory  that  has  yet  to  be  fully  recognized.  Although 
results  are  promising,  these  earth-orbiting  tori  methods  are,  in  a  fashion,  disconnected  from 
the  actual  theory.  This  is  true  because  they  do  not  utilize  the  original  method  used  to  prove 
the  theory.  That  method  uses  a  canonical  transformation  between  two  systems  that  exist  in 
action-angle  coordinates.  So,  the  original  system  must  possess  a  full  range  of  frequencies, 
termed  non-degenerate.  The  typical  approach  in  perturbation  theory  is  to  use  the  two  body 
problem  as  the  original  system,  which  is  fully  degenerate  since  two  of  the  three  fundamental 
frequencies  are  zero.  Thus,  this  does  not  contain  any  range  of  frequencies  and  is  therefore 
not  acceptable.  So,  in  these  techniques,  near  realistic  orbits  were  simulated  or  actual  data 
was  received  and  frequencies  were  extracted  after  the  fact  through  various  spectral  methods 
[23], 

Unlike  previous  attempts,  using  Vinti’s  solution  as  a  starting  point  is  encouraging 
in  that  it  produces  the  full  range  of  frequencies  needed  to  perform  necessary  canonical 
transformations.  The  promise  of  frequencies  is  due  to  the  perturbing  zonal  harmonics 
accounted  for  therein.  Thus,  having  those  frequencies  as  a  result  of  the  solvable  solution 
allows  for  perturbation  theory  to  be  applied  with  respect  to  these  frequencies.  The 
significant  benefit  in  such  a  demonstration  is  the  potential  to  easily  model  long  term 
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behavior  of  the  motion  due  to  highly  accurate  basis  frequencies.  Slight  frequency 
inaccuracies  exist  in  modern  OD  software,  such  as  SGP4,  which  tends  to  a  divergent 
solution  after  a  period  of  time.  However,  if  the  frequencies  can  match  the  dynamics 
solution,  the  current  position  error  will  not  grow  with  time  but  will  slightly  oscillate  around 
the  truth  within  a  reasonable  margin. 

2.9  Summary 

With  the  growth  of  orbital  debris  and  the  expected  significant  increase  of  catalogued 
items,  accurate  analytical  solutions  should  be  sought  after  with  the  fervor  in  which 
they  were  in  the  late  1950s.  Complacency  in  traditional  methods  has  the  potential 
for  catastrophic  results  by  creating  millions  of  conjunction  warnings  when  only  two  or 
three  require  attention.  The  stage  is  set  to  explore  other  avenues  with  which  to  pursue 
advancement  in  this  field.  What  follows  is  an  examination  of  one  of  these  options.  The 
solution  used  is  only  at  the  beginning  of  its  development  but  will  be  demonstrated  to 
already  compete  with  SGP4  in  a  perturbation-limited  comparison.  A  method  to  begin 
adding  perturbations  to  this  solution  is  explored  by  accounting  for  the  next  largest  effect 
past  the  geopotential,  air  drag. 


28 


III.  Orbit  Fitting 


This  chapter  details  how  the  current  research  will  leverage  the  newly  revised  Vinti 
solution  in  an  orbit  fitting  scheme.  A  nonlinear  least  squares  algorithm  is  implemented 
for  differential  correction  of  a  reference  or  estimated  orbit  based  on  provided  observations. 
This  method  uses  the  Vinti  state  transition  matrix  that  is  found  through  finite  differencing 
techniques.  Validation  of  the  matrix  calculation  and  initial  fitting  performance  is  provided. 

3.1  Using  Nonlinear  Least  Squares 

Among  the  various  estimation  techniques  available,  the  current  research  will  utilize 
the  capabilities  of  a  batch  method,  nonlinear  least  squares  routine  to  perform  orbit  fitting 
using  a  set  or  “batch”  of  observations.  Benefits  inherent  in  least  squares  that  are  attractive 
for  the  current  research  include  simplicity,  stability,  and  speed.  On  the  other  hand,  Kalman 
filters  would  be  best  suited  for  sequential  estimation  of  stochastic  or  random  components 
of  a  dynamics  or  measurement  model  where  the  prediction  is  improved  upon  continuously 
with  each  piece  of  new  data. 

This  research  hopes  to  minimize  the  amount  of  random  behavior  and  drill  down 
to  the  core  of  a  Vinti  model  examination.  When  considering  orbital  motion,  the  main 
perturbing  effects  that  are  the  least  predictable  include  those  of  atmospheric  fluctuations 
(to  be  discussed  in  Chapter  V)  due  to  solar  activity  and  the  like.  Further,  a  batch  method 
is  used  to  mimic  how  JSPOC  performs  orbit  fitting.  Observations  worldwide  for  a  specific 
catalog  item  are  processed  each  time  a  specific,  analytical  orbital  solution  is  improved 
upon. 

A  more  rigorous  derivation  and  statistical  treatment  of  this  method  can  be  found  in 
Vallado  [27:728-753]  and  Wiesel  [75:60-73].  What  follows  describes  the  relevant  aspects 
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of  the  overall  algorithm.  Key  summaries  are  included  here  for  clarity  and  to  identify 
necessary  quantities. 

3.1.1  An  Introduction 

We  begin  with  a  system  whose  dynamics  are  shown  as  equations  of  motion 

X(t)  =  f(X,  t)  (3.1) 

to  be  numerically  integrated  or  as  an  explicit  solution 

X(f)  =  h(X(f0),  0-  (3.2) 

where  X  =  (x,y,z,  x,y,z)T.  This  nonlinear  mapping  in  Equation  3.2  best  represents  the 
Vinti  solution:  one  directly  obtains  the  state  X  at  some  time  t,  given  the  state  at  time  to- 
The  dynamics  of  the  orbital  motion  are  assumed  to  be  deterministic  without  random  noise. 
However,  as  is  the  fundamental  requirement  for  all  estimation  techniques,  the  true  states 
are  unknown  and  must  be  estimated. 

In  this  method,  a  reference  trajectory  is  required.  State  estimates  calculated  from  this 
guess  of  an  orbit  serve  as  predictions  for  where  we  believe  the  satellite  will  be  at  some  time 
in  the  future.  The  goal  of  our  estimation  scheme  is  to  modify  this  reference  orbit  to  be  as 
close  to  reality  as  possible. 

To  begin,  a  state  vector  at  some  epoch  time  is  needed  to  describe  this  orbit.  In  catalog 
maintenance,  this  reference  trajectory  results  from  a  previous  orbit  fit,  perhaps  one  day 
or  one  week  in  the  past.  The  current  batch  of  data  to  process  represents  all  observations 
received  since  the  last  fit.  Characterizing  maneuvers  of  active  satellites  that  have  occurred 
since  the  last  fit  or  other  external  circumstances  is  outside  the  of  scope  of  the  current 
research. 

From  the  perspective  of  catalog  maintenance,  the  state  at  epoch  can  be  taken  from  a 
previous  fit.  For  an  uncorrelated  track,  or  observations  from  an  orbit  that  does  not  readily 
line  up  with  another  cataloged  item,  this  reference  guess  is  generated  using  some  initial 
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determination  scheme.  This  can  be  a  complex  process  and  will  not  be  covered  here.  It  is 
then  corrected  to  result  in  a  trajectory  that  is  satisfactorily  close  to  the  truth. 

Deterministic  dynamics  allows  for  the  expression 

6X(t)  =  Q>(t,to)6X(t0)  (3.3) 

to  be  valid,  where  6X  is  some  small  variation  to  the  state.  <1>(7,  to)  is  termed  the  state 
transition  matrix  and  is  a  Jacobian  matrix  containing  partials  derivatives  of  the  state  X  at 
time  t  with  respect  to  the  state  at  time  to,  to  be  found  later  (see  §3. 1.5.1).  This  allows  us  to 
say  that  some  small  displacement  at  time  /0  can  be  propagated  to  some  time  t.  Effectively, 
this  shows  us  how  the  state  at  time  t  will  be  affected  by  a  minor  adjustment  at  the  epoch 
time.  This  concept  becomes  important  when  we  try  to  minimize  the  difference  between  our 
predicted  orbit  and  the  one  being  observed. 

A  numerically  integrated  orbit  serving  as  a  truth  model  (see  §3.1.3)  generates  a  batch 
of  observations,  z,  to  be  used  in  fitting.  This  z  vector  is  usually  related  to  the  state  vector 
through  some  nonlinear  function  G(X).  To  simplify  analysis,  this  relationship  will  be  one 
to  one  with  the  position  vector  of  the  truth  orbit  or  z  =  G(X)  =  (x,y,z)Jruth.  Therefore,  no 
complex  linearizations  or  observation  relationships  are  needed.  Also,  weights  that  account 
for  noise  expected  in  the  measurement  will  not  be  used  but  could  easily  be  incorporated  for 
realistic  applications.  The  general  process  follows. 

3.1.2  The  Algorithm 

Derivation  of  the  linear  form  of  least  squares  gives  the  normal  equation 

Xffo)  =  ( TtQTx  T)~l  TtQ~ 1  z  (3.4) 

where  one  can  directly  solve  for  the  state  given  a  set  of  observations.  Note  that  the  matrices 
on  the  right  are  accumulated  values  and  the  z  vector  contains  all  measurements  from  the 
batch.  However,  with  the  high  nonlinearity  of  the  problem,  it  is  necessary  to  solve  for 
corrections  to  a  predicted  orbit  based  on  residuals.  This  equation  is  identical  in  form  after 
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these  quantities  are  exchanged 


SX(to)  =  (: TtQT1  T)~l  TtOT{  r 


(3.5) 


Above,  Q  represents  a  weighting  matrix  that  can  account  for  noisy  measurements.  T  is  the 
linearization  of  the  measurement  multiplied  by  the  state  transition  matrix  and  can  be  found 
as 


T  =  HdXtJo) 


(3.6) 


where 


1  0  0  0  0  0 
0  1  0  0  0  0 
0  0  1  0  0  0 


(3.7) 


Since  no  process  noise  will  be  added  to  the  measurements,  Q  will  be  taken  as  the  identity 
matrix.  However,  it  will  remain  in  the  explanation  below  for  generality. 

A  step-by-step  explanation  of  the  iterative  process  will  now  be  provided  as  well  as 
a  flow  diagram,  see  Figure  3.1.  Given  an  initial  guess  of  the  epoch  position  and  velocity 
defining  the  reference  orbit,  perform  the  following  for  each  observation  time,  tj,  for  all  n 
observations: 


•  Using  the  Vinti  package7,  produce  a  state  estimate  and  calculate  the  state  transition 
matrix,  0(f,-,  to) 

•  Subtract  the  predicted  position  from  the  observed  state  producing  the  residual  vector 

A  A  —  HXjprejiCte(i 

•  Accumulate  values  to  running  sums  of  the  matrix 

i 

YjTJQ-'T,  (3.8) 

i=0 

7  Vinti  package  refers  to  the  set  of  routines  that  calculates  Vinti’s  solution  as  modified  by  Wiesel.  It  inputs 
a  state  vector  at  some  epoch  time  and  outputs  various  orbital  parameters  at  some  other  specified  time. 
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and  the  vector 


i= 0 


After  all  n  observations  are  processed  (when  i  =  n). 


(3.9) 


The  correction  to  the  reference  orbit  can  be  calculated  as 

(  n  A  1  n 


dX(to)  - 


V  i=0 


YjT?Qi'  Ti  YjTjQ-'r, 


) 


i= 0 


(3.10) 


Correct  the  reference  orbit 


Xre/'(?i o)  -  Xre/(to)  +  dX(?0) 


(3.11) 


If  convergence  criteria  is  met,  the  process  is  complete  with  Xref  being  our  best 
estimate  of  the  true  orbit.  If  not,  begin  again  to  successively  calculate  corrections. 
Convergence  criteria  is  unique  to  each  application  but  is  essentially  a  way  to 
determine  when  the  estimation  routine  has  exhausted  its  ability  to  find  a  better 
solution.  Using  a  parameter  to  indicate  goodness  of  a  fit,  Root  Mean  Square  (RMS, 
to  be  explained  shortly)  can  be  compared  between  iterations  to  see  if  the  solution  is 
improving  and  if  the  routine  should  continue.  Determining  the  percentage  change, 
criteria  to  quit  can  be  set  with 


RMS  „id  -  RMS  new 
RMS  0id 


(3.12) 


Convergence  during  this  research  was  typically  declared  after  only  3  or  4  iterations 
while  maximum  iterations  were  set  to  10.  Setting  a  maximum  number  is  a  second 
way  to  decide  when  to  quit  by  admitting  that  the  solution  will  not  improve  by  simply 
repeating  the  process  endlessly. 


Examine  and  analyze  residuals. 


Figure  3.1:  Iterative  Nonlinear  Least  Squares  Fitting  Routine 


A  Root  Sum  Square  (RSS)  of  the  residuals  can  be  used  to  plot  and  visualize  model 
divergence,  or  error  growth,  over  time.  With  this,  the  three  components  of  the  residual 
vector  for  each  observation  are  squared,  summed,  then  the  square  root  is  taken.  This  looks 
like 

RSS  j  =  |A/r2  +  r22  +  r33j  =  |r,|  (3.13) 

This  represents  the  magnitude  of  the  estimated  error  at  that  time.  Then,  for  n  observations, 
these  values  can  then  be  averaged  over  the  entire  fit  to  give  an  average  error, 

1  " 

RSSavg  =  -  YW  (3.14) 

nU 

Root  Mean  Square  (RMS)  error  is  the  basis  for  least  squares  methods  due  to  its 
numerical  stability  and  speed.  Residual  components  are  squared  and  summed  for  each 
observation.  When  this  is  complete  for  all  observations,  the  mean  is  taken  and  then  the 
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square  root  as  in 


RMS  = 


A 


1  " 


2  +  r2  +  r2 


), 


(3.15) 


1=0 


This  is  computationally  more  efficient  than  RSS  because  the  square  root  operation  is 
required  only  once  at  the  end  of  data  accumulation  as  opposed  to  once  every  observation. 
RMS  can  provide  a  snapshot  of  fit  performance  although  it  is  not  the  mathematical  average 
error.  In  practice,  RMS  and  RSS  values  are  similar  in  magnitude  and  are  related  by 


RMS  =  -~RSSavg  (3.16) 

yn 

3.1.3  Truth  Model 

Before  moving  on  with  the  development  of  the  fitting  routine,  a  description  of  the 
truth  model  is  in  order.  A  numerical  integration  routine  using  a  fourth-order  predictor 
corrector  algorithm  is  used.  It  uses  a  Hamming  integrator  that  is  known  for  stability 
and  speed.  Flexibility  is  inherent  with  this  tool  as  various  geopotential  scenarios  can 
be  implemented.  For  instance,  a  zonals  only  or  full  geopotential  (to  specified  order  and 
degree)  can  be  declared.  Unique  to  the  current  research,  a  Vinti  geopotential  can  also  be 
invoked  simulating  the  exact  potential  that  is  accounted  for  in  Vinti’s  solution.  This  is 
used  in  initially  validating  the  fitting  scheme  and  later  in  comparing  drag-only  perturbation 
methods.  In  Chapter  V,  the  process  of  using  this  truth  model  to  estimate  air  drag  effects 
will  be  explained. 

3.1.4  Dimensionless  Units 

When  using  orbital  velocities  and  distances,  units  can  be  problematic.  For  example,  if 
metric  units  are  desired  a  low  earth  satellite  may  have  a  radius  vector  on  the  order  of  7000 
km  with  a  velocity  of  6  or  7  km/sec.  Maneuvers  or  changes  in  velocity  would  be  on  the 
order  of  centimeters  or  meters  per  second.  Numerical  conditioning  can  quickly  become  an 
issue  when  dividing  a  large  number  by  a  very  small  number. 
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This  research  will  utilize  geocentric  dimensionless  distance  units  (DU)  and  time  units 
(TU).  These  result  from  scaling  position  and  time  relative  to  the  size  and  rotational  speed 
of  the  earth.  DUs  are  calculated  by  dividing  by  the  average  radius  of  the  earth.  For  metric 
units,  this  means  1  DU  ~  6378.1363  m.  TUs  are  found  by  determining  the  speed  at  which 
an  object  would  orbit  the  center  of  the  earth  from  this  average  radius  but  at  the  equator.  This 
is  found  by  the  expression  ™  and  /j  is  earth’s  gravitational  constant  multiplied  by 

its  mass.  This  works  out  to  be  1  TU  ~  806.81  sec.  Using  these  nondimensionalized  units 
results  in  orbital  parameters  on  the  order  of  unity  which  is  preferred  for  numerical  reasons. 

3.1.5  State  Transition  Matrix 

Obviously,  with  the  estimation  approach  laid  out  above,  finding  the  partials  matrix 
or  state  transition  matrix  using  Vinti’s  dynamics  is  critical.  The  information  provided  by 
the  Phi  matrix  can  be  visualized  as  providing  the  slope  of  the  solution  space  around  the 
reference  orbit  with  which  to  modify  it,  creating  a  more  accurate  orbit  fit.  Considering 
the  numerous  changes  of  variables  and  chain  rule  operations  associated  with  taking  these 
partial  derivatives  analytically,  they  are  performed  numerically.  This  approach  is  called 
finite  differencing  and  is  not  uncommon  especially  with  rapid  analytical  propagators. 
Vallado  says  the  time  and  money  needed  to  develop  elegant  analytical  solutions  is  most 
likely  not  worth  the  slight  overhead  needed  to  simply  go  ahead  and  calculate  them  through 
numerical  partial  differentiations  [76]. 

3. 1.5.1  Finite  Differencing 

The  matrix  of  interest  is  a  Jacobian  and  can  be  displayed  as 


dx(t ) 
dx(t0) 

dx(t) 
dy(t o) 

dx(t ) 
dz(to) 

dX 

®(t,t o)  =  ,)V  = 

OAq 

dy(t) 

dx(t0) 

dy(t) 

dy(t0) 

dy(t) 

dz(to) 

(3.17) 

dz(t) 

dx(t0) 

dz(t) 

dy(t0) 

dz(t) 

<?z0o) 
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First,  the  reference  orbit  is  used  to  generate  several  other  orbits  using  finite  displacements 
in  each  direction  of  the  position  and  velocity  at  the  epoch  time.  This  is  repeated  for  -6Xj, 
+2 6Xj,  and  -28  Xt  for  a  total  of  24  additional  orbits.  These  displaced  orbits  are  then  run 
through  the  Vinti  solution  package  resulting  in  position,  velocity  vectors  at  the  observation 
time.  Figure  3.2  illustrates  this  process  and  Figure  3.3  shows  what  this  would  look  like  in 
configuration  space  for  three  position  displacements. 


Figure  3.2:  Construction  of  Nearby  Orbits 


Figure  3.3:  Visualization  of  Three  Nearby  Orbits  (Configuration  Space) 
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The  elements  of  the  Phi  matrix  are  then  calculated  following 

*  Xx  (t2^X  ~  26'^  ~  \^X  ~  6X)  +  \^X  +  6X)  ~  I \f^x  +  26x))  (3-18) 

a  five  point  Lagrangian  partial  derivative  [77:25.3.6].  This  method  can  be  likened  to  how 
an  Unscented  Kalman  Filter  samples  a  system’s  nonlinearities  directly,  avoiding  linearizing 
the  dynamics  in  the  partials  matrices.  When  nonlinear  effects  become  large  for  certain 
orbits,  an  added  performance  boost  may  be  expected  by  avoiding  linearization  errors. 

The  magnitude  of  the  position  displacements  is  determined  by  taking  the  magnitude  of 
the  position  vector,  in  DU,  and  multiplying  by  10-5.  The  same  process  is  accomplished  for 
velocity  displacements  relative  to  the  velocity  vector,  in  DU/TU.  This  value  was  chosen  so 
it  makes  the  position  and  velocity  state  displacements  on  the  order  of  75  m  and  7.5  cm/sec, 
respectively.  Using  these  dimensionless  parameters,  this  requires  5  or  6  significant  figures 
and  with  double  precision,  there  are  digits  to  spare.  Any  larger  and  there  would  not  be 
sufficient  granularity  required  to  fully  characterize  the  orbit.  Any  smaller  and  there  would 
be  risk  of  round-off  error.  Further,  using  geocentric,  nondimensionalized  units  of  DU  and 
TU  avoids  numerical  problems  encountered  when  dividing  large  astronautical  distances  by 
small  velocities. 

3. 1.5.2  Validation 

To  ensure  the  state  transition  matrix  at  hand  is  usable  and  accurate,  it  must  be 
validated.  Using  linear  systems  theory,  the  Phi  matrix  follows  certain  properties.  First, 
if  the  predicted  time  equals  the  epoch  time,  the  result  is 

O(roUo)  =  /  (3.19) 

Next,  the  change  in  behavior  over  time  should  follow 

6  =  AO  (3.20) 
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where  A  is  a  6x6  matrix  and  when  broken  down  into  the  four  3x3  sub-matrices  looks  like 


0 

I 

^2,1 

0 

(3.21) 


A2, i  contains  the  gravitational  acceleration  terms.  The  upper  right  block  should  be  an 
identity  matrix  due  to  state  space  relation  between  position  and  velocity.  The  diagonal 
blocks  should  be  nulls  as  a  result  of  absence  of  these  terms  in  the  equations  of  motion. 

To  test  this  structure,  an  incremental  predicted  time  supplied,  or  Olbf,  to)  should  result 
in  the  form 


Ofa  +  8t)  -  Ofo) 
St 


Al 


(3.22) 


These  checks  were  performed  adequately  where  the  resulting  matrices  were  of  the 
necessary  structure. 

Next,  an  idea  of  how  long  the  difference  approximation  of  the  Phi  matrix  would  be 
valid  is  examined.  For  this,  a  truth  model  Phi  matrix  was  used  for  comparison.  In  the 
numerical  integration  truth  model  used,  the  capability  exists  to  also  integrate  Equation  3.20 
above.  The  matrices  from  both  models  were  calculated  incrementally  over  a  period  of  time 
and  compared.  At  each  time  step,  the  difference  was  taken  creating  what  can  be  termed  a 
del  matrix.  Various  norms8  were  then  taken  on  this  del  matrix  and  plotted  over  time.  This 
was  performed  for  a  certain  size  and  inclination  of  orbit9  and  then  repeated  over  various 
eccentricities.  Results  from  a  few  of  these  are  shown  in  Figures  3.4  to  3.6. 

The  Phi  matrix  as  found  through  finite  differentiation  seems  to  remain  within  a 
reasonable  margin  of  the  integrated  one  over  a  two  day  span.  As  illustrated,  some  errors 
oscillate  within  the  same  tolerance  and  others  grow  apart  over  time.  Either  way,  the  value 
of  the  difference  is  less  than  6%  over  the  time  examined  except  for  when  eccentricity 
approaches  0.9.  Considering  the  time  periods  of  interest  in  orbit  fitting  and  the  purposes 

8 1-norm:  Max  column  sum;  Frobenius  norm:  Square  root  of  the  sum  of  the  square  of  all  elements; 
Infinity-norm:  Max  row  sum 

9altitude  at  perigee:  ~  3,500  km,  inclination:  20° 
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of  this  research,  the  Phi  calculations  are  implemented  appropriately  to  be  inserted  into 
the  fitting  routine.  The  Phi  quantities  will  still  provide  the  slope  in  solution  space  for 
convergence,  however  quadratic  convergence,  common  in  least  squares  estimation,  may  be 
lost  and  require  more  iterations  to  complete  the  routine. 


Figure  3.4:  Norms  of  Phi  Difference:  e  =  0.001 


3.2  Initial  Orbit  Fitting  Performance 

With  the  Phi  matrix  available  and  validated,  the  least  squares  orbit  fitting  scheme 
can  be  tested.  As  an  initial  check  of  the  model  and  code,  sample  observations  from  a 
trajectory  perturbed  only  by  Vinti’s  geopotential10  are  created  every  minute  for  one  day  and 
supplied  to  the  Vinti  OD  software.  This  limited  perturbation  approach  allows  for  an  initial 
comparison  between  truth  and  estimation  methods  by  keeping  the  observation  perturbations 
equal  to  those  accounted  for  in  Vinti’s  solution.  The  only  differences  should  be  numerical 
and  can  help  characterize  errors  in  the  model  or  calculations. 

10zonal  potential  terms  past  Ji  are  approximated  as  J4  -  -J\,  J(,  =  J[  and  so  on 
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Figure  3.5:  Norms  of  Phi  Difference:  e  =  0.1 


Figure  3.6:  Norms  of  Phi  Difference:  e  =  0.7 


Five  sample  orbits  with  various  parameters  were  created  and  fit  using  the  above 
methodology.  Results  are  summarized  in  Table  3.1.  As  expected,  performance  varies 
depending  on  size  and  shape  of  the  orbit.  However,  centimeter-level  fits  provide  confidence 
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that  the  dynamics  model,  Phi  calculations,  and  fitting  mechanisms  are  working  in  concert 
as  expected.  Wiesel’s  advertised  error  is  commensurate  with  these  magnitudes  [23]. 


Table  3.1:  Vinti  Geopotential  Orbit  Fits 


Test  Orbit 

Perigee  Altitude  (km) 

e 

i 

RMS  (cm) 

1 

400 

0.01 

28.5° 

29.5 

2 

500 

0.2 

45° 

20.8 

3 

800 

0.2 

28.5° 

17.9 

4 

1000 

0.7 

28.5° 

74.6 

5 

1000 

0.001 

0.01° 

9.8 

To  provide  more  insight  into  the  RMS  values,  residuals  over  time  are  shown  for  the 
first  two  orbits  in  Figures  3.7  and  3.8.  These  illustrations  decompose  where  the  errors 
manifest,  allowing  understanding  into  how  the  model  fails  to  match  reality.  Again,  as 
expected,  these  differences  seem  to  be  a  function  of  the  size  and  shape  of  the  orbit.  There 
appears  to  be  a  bias  in  the  radial  direction  for  all  orbits  tested.  This  bias  is  speculated 
to  result  from  a  period  compensation.  Adjusting  radial  position  by  the  model  is  the  most 
direct  way  to  make  up  for  inaccuracies  in  modeling  the  orbital  period.  However,  for  Orbit 
1  (Figure  3.7),  there  appears  to  be  a  more  pronounced  difficulty  in  matching  the  period 
as  evidenced  by  the  linear  growth  on  either  side  of  the  midpoint  of  the  radial  and  in¬ 
track  residuals.  Otherwise,  the  residuals  seem  to  exhibit  stable  behavior  and  can  possibly 
be  explained  by  noise  in  the  calculations.  Noise  on  the  centimeter  level  is  acceptable 
considering  it  is  on  the  order  of  solution  error.  Three  other  orbits  are  analyzed  and  similar 
figures  can  be  found  in  the  appendix  for  reference. 
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Figure  3.7:  Vinti  Geopotential  Fit  Residuals:  Orbit  1 


Figure  3.8:  Vinti  Geopotential  Fit  Residuals:  Orbit  2 


3.3  Conclusion 

This  chapter  demonstrated  the  orbit  fitting  method  for  the  current  research  to 
include  finding  the  state  transition  matrix.  Preliminary  checks  provide  confidence  in  the 
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implementation  and  development  of  computer  code  to  perform  needed  calculations  for 
further  analysis.  Initial  performance  of  the  fully  implemented  fitting  procedures  lays  the 
foundation  upon  which  to  examine  how  well  the  Vinti  solution  can  approximate  the  full 
geopotential  and  effects  due  to  air  drag. 
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IV.  Results 


This  chapter  reviews  the  performance  of  the  Vinti  fitting  routine  as  developed  in 
Chapter  III  when  used  to  fit  orbits  perturbed  by  the  full  geopotential.  A  wide  survey  of 
fit  performance  is  taken  across  various  inclinations  and  eccentricities.  Next,  a  detailed 
analysis  is  given  for  two  representative  orbits.  These  demonstrations  help  identify  best 
case  and  worst  case  behavior  within  a  functional  orbital  regime  for  this  solution.  These 
same  example  orbits  are  then  used  for  an  exercise  in  predicting  future  miss  distances. 

4.1  Fitting  Performance 

Following  the  procedure  explained  in  Chapter  III,  truth  observations  are  created  once 
every  minute  for  one  day  using  different  sizes  and  shapes  of  orbits.  Unlike  earlier,  these 
orbits  are  perturbed  by  a  full,  20x20  geopotential.  For  added  perspective,  these  orbits  are 
also  provided  to  a  least  squares  fitting  routine  using  an  unclassified  version  of  SGP4.  Using 
the  same  test  orbits  from  Table  3.1,  Table  4.1  shows  the  resulting  RMS  values.  Again,  RMS 
indicates  goodness  of  a  fit. 

Table  4.1:  RMS  Values  For  Full  Geopotential  Orbit  Fits  (m) 


Test  Orbit 

Vinti 

SGP4 

1 

603 

500 

2 

220 

878 

3 

245 

982 

4 

108 

5154 

5 

458 

458 
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At  first  glance  when  comparing  RMS  values,  Vinti  seems  to  outperform  SGP4, 
however,  it  should  be  noted  that  SGP4  is  at  a  slight  disadvantage.  The  truth  observations 
are  only  perturbed  by  a  full  geopotential,  while  SGP4  is  attempting  to  incorporate  all 
perturbations  with  the  exception  of  air  drag.  This  disconnect  becomes  more  prominent 
when  the  apogee  of  an  eccentric  orbit  is  high  enough  for  third  body  effects  to  be  noticeable, 
Orbit  4  for  example  which  reaches  an  altitude  of  over  35,000  km.  Therefore,  this  does 
not  provide  an  ultimate  performance  comparison  between  the  two  models.  What  it  does 
provide,  though,  is  some  relative  metric  from  which  to  draw  some  inference.  In  the 
future,  these  additional  perturbations  could  be  added  to  the  solution  for  a  more  accurate 
representation  of  reality. 

Even  considering  the  disadvantage  of  a  slight  mismatch  of  perturbations,  the  fact  that 
the  Vinti  model  is  mostly  doing  better  than  SGP4  is  promising.  The  geopotential  is  the 
dominant  perturbation  in  these  low  earth  orbits.  SGP4  is  at  an  end  when  approximating 
for  these  effects  by  truncating  after  expansions  through  J5.  However,  the  Vinti  model  is 
only  using  the  solvable  solution  which  accounts  exactly  for  J2  and  nearly  three-fourths  of 
/4.  The  true  test  will  be  when  the  remaining  geopotential  effect  is  incorporated  through  a 
general  perturbations  approach.  This  is  left  for  future  works. 

4.1.1  A  Wide  Survey  of  Orbits 

Before  detailed  information  of  any  one  orbit  is  presented,  a  wide  range  of  orbits  is 
sampled.  This  gives  an  idea  of  relative  best  case/worst  case  scenarios  while  illustrating 
where  the  model  works  and  where  it  struggles.  Figure  4.1  shows  a  sampling  of 
eccentricities  and  inclinations  for  a  perigee  height  of  400  km. 

It  is  evident  the  Vinti  solution  as  modified  by  Wiesel  begins  to  break  down  for 
inclinations  greater  than  60°  or  so.  However,  this  is  not  as  a  result  of  any  critical  inclination 
difficulties  in  the  sense  of  singularities  resulting  from  a  zero  apsidal  rate.  All  data  indicates 
singularities  do  not  exist  at  or  near  63.4  as  there  does  in  other  theories  such  as  Brouwer’s 
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original  solution.  These  errors  seem  to  grow  to  infinity  when  approaching  the  polar  regime. 
When  increasing  inclination,  an  RMS  of  around  25  km  was  calculated  at  82.5°.  When 
decreasing  inclination  towards  polar,  similar  RMS  was  found  at  97.5°.  A  solution  does  not 
exist  between  due  to  a  numerical  flaw  in  the  Vinti  integrals.  Resolving  this  flaw  is  outside 
the  scope  of  the  current  research. 
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Figure  4.1:  Orbit  Fit  Performance:  Perigee  Height  =  400  km 


When  inclination  is  increased  past  polar,  the  errors  eventually  come  down  and 
behavior  begins  to  resemble  what  it  was  at  inclinations  below  the  first  critical  inclination. 
This  occurs  near  the  other  critical  inclination  which  is  about  116.6°  (for  reference,  the 
critical  inclinations  are  noted  on  Figure  4.2  with  blue  arrows  and  the  reader  is  directed 
to  Appendix  C  for  more  details  on  the  critical  inclination  phenomena).  The  argument  of 
perigee  is  fixed  in  inertial  space  at  these  two  locations  but  apsidal  rate  is  negative  between 
them.  It  is  the  author’s  belief  that  this  negative  rate  is  the  source  of  the  numerical  flaw 
within  the  Vinti  integrals  as  modified  by  Wiesel. 

Vinti’s  original  solution  was  claimed  to  be  accurate  through  high  inclination.  This 
means  a  negative  rate  should  not  have  been  an  issue.  The  difference  is  suspected  to  manifest 
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in  Wiesel’s  root  finding  scheme  illustrated  in  §2.7.  To  solve  for  some  of  the  terms  in 
his  elliptic  integrals  numerically,  he  integrates  from  perigee  to  apogee.  Therefore,  if  this 
negative  rotation  is  not  accounted  for  properly,  a  less  accurate  solution,  and  ultimately  a 
singularity,  could  be  the  result. 

Singularities  do  exist  for  equatorial  and  near  polar  orbits  as  well  as  eccentricity  of 
zero.  Vinti’s  original  theory  did  not  have  such  problems  but  Wiesel’s  numerical  solving 
routine  relies  on  certain  root  finding  methods  that  need  a  definite  perigee  and  northernmost 
point.  However,  eccentricities  as  low  as  0.001  and  inclinations  of  0.01  have  resulted  in 
satisfactory  fits.  Detailed  analysis  for  performance  of  this  fitting  routine  will  focus  on  a 
functioning  orbital  regime  of  non-circular  orbits  in  the  inclination  range  of  greater  than  0° 
but  less  than  65°. 


Figure  4.2:  High  Inclination  Issue:  Perigee  Height  =  400  km,  Eccentricity  =  0.2 

4.1.2  A  Closer  Look 

When  comparing  data  from  Table  4.1  and  Figure  4.1,  there  are  noticeable  variations 
in  performance  for  different  orbital  parameters.  The  first  two  test  orbits  in  the  table  are 
identified  as  candidates  for  further  analysis.  They  are  used  because  they  seem  to  represent 
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a  worst  case  (Orbit  1)  and  best  case  (Orbit  2)  scenario  from  the  functional  range  of 
inclinations  with  a  reasonable  eccentricity  of  0.2  or  lower. 

Figure  4.3  displays  residuals  for  the  one  day  fit.  Since  the  observations  are  created 
from  the  truth  model  with  no  process  noise,  the  residuals  are  actually  the  error.  These 
terms  will  be  used  interchangeably. 


Figure  4.3:  Full  Geopotential  Fit  Residuals:  Orbit  1 


Notice  all  the  residuals  exhibit  a  periodic  behavior  that  coincides  with  the  orbital 
period.  In  addition,  the  in-track  error  has  an  extra,  longer  oscillation  with  a  period  of 
about  half  of  a  day.  For  orbit  1  which  has  an  RMS  of  603  m,  Table  4.2  shows  more  detailed 
statistics  about  this  fit.  Considering  the  scale,  all  the  errors  are  near  zero  mean  so  there 
is  no  bias  within  the  day  fit.  When  the  data  is  fit  to  a  normal  distribution,  the  standard 
deviation,  cr,  illustrates  how  spread  out  the  errors  are.  This  combined  with  the  minimum 
and  maximum  errors  show  the  in-track  errors  are  by  far  dominant.  This  artifact  was  also 
illustrated  by  Gordon  in  1978  when  he  compared  Vinti,  Brouwer,  and  Brouwer- Lyddane 
solutions  [52]. 
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Table  4.2:  Error  Statistics  From  Orbit  1  Fit  (m) 


Average 

cr 

Min 

Max 

Radial 

-4.6 

180 

-467 

411 

In-track 

0.03 

554 

-1241 

1334 

Orbit  Normal 

-5.9 

152 

-355 

347 

Orbit  2  is  a  much  better  fit  with  an  RMS  value  of  220  m.  Figure  4.4  shows  a  similar 
periodic  trend  but  the  longer  frequency  error  is  less  pronounced  in-track.  Again,  Table  4.3 
shows  that  means  were  near  zero  and  the  error  in  the  velocity  direction  was  the  greatest 
contributor  to  overall  error. 


Figure  4.4:  Full  Geopotential  Fit  Residuals:  Orbit  2 


4.2  Error  Growth 

Wiesel’s  solution  at  the  heart  of  a  fitting  routine  has  displayed  initial  promise  for 
adequate  results  for  a  certain  range  of  orbital  parameters.  It  is,  at  the  least,  comparable 
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Table  4.3:  Error  Statistics  From  Orbit  2  Fit  (m) 


Average 

<T 

Min 

Max 

Radial 

-1.3 

50 

-147 

126 

In-track 

-0.2 

201 

-353 

456 

Orbit  normal 

2.6 

76 

-181 

141 

to  how  SGP4  performs  when  fitting  the  same  data.  However,  fitted  data  is  in  the  past  and 
potential  collisions  are  in  the  future.  Therefore,  an  orbital  propagator  is  only  as  good  as 
how  well  it  can  fit  an  orbit  then  predict  a  future  position/velocity  state.  This,  combined 
with  other  fits  and  predictions,  provide  miss  distances  for  nearby  orbiting  craft  or  debris. 

An  idea  of  how  long  a  solution  will  be  valid  is  explored.  Data  is  fit  for  one  day  then 
used  to  predict  future  positions.  This  is  compared  to  what  the  integrated  truth  model  outputs 
for  those  future  times  using  the  same  orbit.  Orbits  1  and  2  from  above  are  examined  again 
in  this  exercise. 

After  the  magnitude  of  the  position  error  at  each  timestep  is  calculated,  Figures  4.5  and 
4.6  show  the  behavior  in  future  predictions  as  compared  to  the  integrated  truth.  Notice  the 
difference  in  scale  between  the  two.  In  both  cases,  an  increasing  amplitude  periodic  error 
exists  while  a  secular  growth  is  evident.  When  considering  geopotential  effects,  even  zonal 
terms  contribute  to  periodic  changes  in  semi-major  axis  and  eccentricity,  secular  growth  in 
the  node  and  argument  of  perigee,  and  a  combined  secular/periodic  effect  on  mean  anomaly. 
However,  sectoral  and  tesseral  contributions  are  only  periodic  for  all  the  elements.  The 
approximation  of  even  zonals  beyond  J2  and  lack  of  the  remaining  sectoral/tesseral  terms 
can  be  implicated  by  this  error  growth  behavior. 

When  a  linear  fit  is  calculated  for  the  future  error,  the  slope  of  that  line  is  a  good 
estimate  of  the  average  error  growth  rate.  For  the  worst  case,  Orbit  1,  the  average  growth 
is  a  little  over  1  km  per  day.  For  the  best  case,  Orbit  2,  the  error  growth  is  estimated  to  be 
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about  370  m  per  day.  This  difference  can  lend  insight  into  how  well  a  day’s  fit  will  translate 
into  future  error  growth.  Since  these  two  sample  orbits  seem  to  bound  the  performance 
within  usable  orbit  regimes,  it  could  be  inferred  that  average  error  growth  for  other  orbits 
should  lie  between  these  two  values. 

However,  while  the  amplitude  of  the  oscillations  increase,  so  does  the  uncertainty 
for  any  miss  distance  estimate.  For  example  with  Orbit  1,  after  2  days  the  average  error 
could  be  calculated  to  be  approximately  3  km.  However,  it  could  be  1.8  km  or  as  bad  as 
5  km.  After  4  days,  the  average  error  is  just  over  5  km  with  a  minimum  of  3  km  and 
maximum  approaching  9  km.  Considering  the  suspected  sources  of  this  error  growth, 
incorporating  remaining  higher  order  geopotential  effects  into  the  solution  should  pay 
dividends  in  increased  solution  accuracy  over  time. 
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Figure  4.5:  Magnitude  of  Error  Growth:  Orbit  1 
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Figure  4.6:  Magnitude  of  Error  Growth:  Orbit  2 


4.3  Processing  Time 

A  crucial  attribute  for  general  perturbations  solutions  such  as  Vinti  is  undoubtedly 
minimal  overhead  needed  for  computer  processing.  The  goal  of  using  Vinti  in  place 
of  numerical  integrators  is  similar  accuracy  with  rapid  processing.  Also,  as  mentioned 
above,  using  numerical  partials  over  analytical  ones  typically  does  not  result  in  noticeable 
penalties.  Although,  in  practice,  this  is  not  the  case. 

Wiesel’s  solution  is  still  an  analytical  one  in  that  it  does  not  have  to  integrate  the 
equations  of  motion  for  each  time  step  along  the  way.  It  does,  though,  perform  an 
integration  to  generate  values  needed  for  a  specific  orbit.  Doing  this  once  then  propagating 
the  state  forward  in  time  is  still  rapid.  However,  the  fitting  routine  utilized  for  this 
research  requires  several  orbital  solutions  to  be  calculated  for  each  time  step  in  order  to 
create  the  numerical  state  transition  for  the  purposes  of  fitting.  This  integration  within 
Wiesel’s  solution  is  calculated  dozens  of  times  for  each  observation.  This  overhead  makes 
it  considerably  slower  than  SGP4  but  on  the  same  order  of  time  it  takes  to  numerically 
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integrate  a  trajectory  of  similar  duration.  Performing  orbit  fitting  using  an  analytical  Phi 
matrix  should  alleviate  this  penalty. 

4.4  Conclusion 

Using  only  the  solvable  part  of  Vinti’s  solution  as  modified  by  Wiesel,  this  chapter  fit 
orbits  perturbed  by  the  full  geopotential.  Performance,  as  compared  to  SGP4,  is  initially 
promising  considering  no  attempts  have  been  made  to  incorporate  perturbations  beyond 
the  Vinti  geopotential  in  an  effort  to  close  the  gap  towards  the  full  geopotential.  However, 
limitations  were  identified  with  high  inclination  and  in  terms  of  speed.  Solutions  near  high 
inclination  were  not  expected  to  produce  poor  results  due  to  the  original  Vinti  solution 
being  advertised  as  free  of  singularities.  Considering  the  popularity  of  such  orbits  in 
practice,  this  is  a  major  drawback  for  operational  application.  With  regards  to  speed,  a 
rapid  analytical  solution  was  sought.  However,  using  numerical  partials  for  orbit  fitting 
slows  the  process  down  enough  where  processing  time  for  fits  are  about  the  same  as 
integration.  This  provides  motivation  to  repeat  fitting  development  with  aide  of  analytical 
partial  matrices  which  is  suggested  for  future  works.  The  next  chapter  explores  a  method 
to  add  perturbations  to  the  Vinti  solution. 
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V.  Air  Drag 


This  chapter  lays  out  an  approach  to  account  for  and  incorporate  air  drag  effects  in 
Vinti  OD.  The  great  advantage  of  possessing  the  closed  form  action-angle  solution  to  the 
system  will  be  exploited  to  perform  classical  perturbation  theory.  Being  able  to  express 
the  dynamics  in  such  a  unique  way  provides  valuable  and  straightforward  manipulation  of 
coupling  effects  when  adding  perturbations. 

An  extensive  review  of  the  literature  reveals  this  approach  has  not  been  applied  to  the 
original  solution  of  Vinti  or  the  modification  due  to  Wiesel.  The  methodology  is  similar 
in  theory  to  what  Wiesel  lays  out  in  his  text  [35:202-205]  and  demonstrates  in  his  analysis 
of  periodic  orbits  with  low  eccentricity  [73:3,10-19].  Small  displacements  to  a  multiply 
periodic  system  are  considered  and  then  a  small  perturbing  acceleration  is  added.  This 
forcing  term  is  then  Fourier-analyzed  and  integrated,  providing  a  function  of  time  for  the 
perturbations.  An  analytical  solution  is  then  available  to  calculate  displacements  due  to 
drag  at  any  given  time. 

5.1  Air  Drag  Considered 

Up  until  now,  our  dynamics  have  been  assumed  to  be  deterministic.  This  is  to  say  that 
there  is  no  random  component  involved,  whether  or  not  we  can  accurately  account  for  it. 
However,  there  is  one  aspect  of  orbital  motion  that  can  be  most  likened  to  having  random 
behavior  and  not  included  yet:  air  drag.  Intuitively  speaking,  it  does  not  seem  important  to 
consider  the  effects  of  air  flowing  over  a  spacecraft  body.  However,  in  space,  air  interacts 
by  amassing  several  small  collisions,  one  particle  at  a  time,  which  is  contrary  to  the  idea 
of  it  “flowing.”  At  such  speeds,  even  the  small  mass  of  air  molecules  impart  a  significant 
perturbing  acceleration  over  time  through  these  tremendous  collisions.  Of  course,  this 
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effect  diminishes  as  orbital  altitudes  increase  and  atmospheric  density  decreases.  For  low 
earth  satellites,  this  is  the  dominant  perturbation  after  the  earth’s  geopotential. 

Modern  efforts  have  been  made  to  develop  somewhat  deterministic  models  for  the 
earth’s  atmosphere.  This  could  provide  information  to  the  dynamics  model  in  planning  for 
and  accurately  quantifying  the  effect  of  air  drag,  ffowever,  the  behavior  of  the  atmosphere 
is  exceedingly  complex  and  dependent  upon  a  number  of  phenomena  that  are  effectively 
unpredictable.  Needless  to  say,  this  aspect  of  our  dynamics  model  contains  the  most  amount 
of  stochastic,  or  random,  behavior.  An  adequate  consideration  of  air  drag  is  needed  to  make 
Vinti’s  solution  relevant  to  modern  OD. 

Several  techniques  have  been  used  through  the  years  to  include  the  dissipative  effects 
of  air  drag  in  orbital  motion  solutions.  The  original  SGP  used  a  simplistic  mean  motion 
approximation  while  SGP4  uses  power  density  functions  [9].  In  1973,  Vinti  published  a 
method  to  include  perturbing  forces  in  his  original  solution  [43].  He  proposed  a  set  of 
Gaussian  variational  equations  to  solve  for  osculating  elements  in  terms  of  the  variables  he 
used.  This  method  examines  changes  to  the  separation  coordinates  and  momenta,  a,  and 
yS Although  similar  in  principal,  the  current  effort  will  work  in  actions  and  angles,  X,  and 
Of.  These  variables  are  unique  to  Wiesel’s  modification  of  Vinti’s  solution  and,  thus,  have 
not  been  examined  previously  as  a  means  to  solve  for  air  drag. 

5.2  Using  Perturbation  Theory 

Demonstration  of  the  underlying  theory  and  solution  derivation  in  this  section  will  be 
displayed  in  the  Cartesian  reference  frame.  This  is  performed  for  clarity  even  though  the 
system  is  solved  through  various  changes  of  variables. 

Recalling  §3.1,  we  can  express  the  equations  of  motion  of  the  unperturbed  Vinti 
system  as 

X  =  f(X,  t)  (3.1) 
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or  as  an  explicit  solution 


X(t)  -  h(X(?o),  t). 


(3.2) 


where  X  =  (x,y,z,  x,y,z)T.  These  relations  of  a  reference  orbit  are  calculated  by  the  Vinti 
solution  as  modified  by  Wiesel.  Now,  write  Hamilton’s  equations  as 


x  =  z 


m 

ox 


(5.1) 


where 


Z 


0  / 
-/  0 


v  ) 

is  the  symplectic  group  matrix  and  preserves  the  properties 


(5.2) 


Z_1  =  Zr  =  -Z 


(5.3) 


Now,  linearize  the  Hamiltonian  around  the  reference  orbit  and  add  a  forcing, 
perturbing  acceleration.  We  can  now  consider  small  changes  to  the  reference  system  as 


SX  =  Z 


d2(H 

~dX? 


dX  +  X, 


pert 


(5.4) 


Including  this  acceleration  in  the  variational  equation  follows  the  general  assumption  in 
perturbation  theory  that  the  perturbation  will  remain  small.  In  the  case  of  air  drag,  this 
assumption  only  fails  during  the  final  stages  of  orbital  decay  when  the  flight  profile  is 
drastically  altered  by  the  exponentially  growing  atmospheric  density.  The  drag  perturbation 
can  be  inserted  here  and  shown  as 


/  \ 

(  \ 

Xdrag  — 

I*drag 
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HZ 

^drag 

V  7 

&drag 
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where 
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1  CrfApref  p 


tag 


W 


2  nz  /  Vcf 

Combine  the  leading  constant  terms  into  a  modified  ballistic  coefficient, 

CrfApre  f 


B*  = 


2m 


(5.6) 


(5.7) 
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where  Cd,  A,  and  m  are  the  satellite’s  drag  coefficient,  frontal  area,  and  mass,  respectively. 
Further,  pref  is  the  atmospheric  density  at  perigee  while  p  is  the  density  at  the  satellite’s 
location.  (See  §5.8  for  details  concerning  the  atmosphere  model)  The  drag  acceleration  is 
now 

adrag  =  -B*— |v|v  (5.8) 

Pref 

This  expression  for  the  drag  acceleration  can  be  found  in  Vallado[27 : 1 14]  and  Wiesel[73:3] . 
For  special  perturbations,  this  set  of  equations  would  be  given  to  numerical  integrators  for 
solving;  that,  however,  would  defeat  the  purpose  of  the  current  research  of  seeking  out  a 
rapid  analytical  solution.  However,  certain  numerical  techniques  can  be  utilized  to  exploit 
their  benefit  while  employing  an  analytical  method.  The  next  section  demonstrates  this. 


5.3  General  Perturbations  in  Action-Angle  Form 

Consider  a  change  of  variables  to  action-angle  form  where  Y  =  (61,62,63, 1  i,I2,I3)t 
and  %  is  the  Hamiltonian  in  terms  of  these  new  variables,  or  more  concisely, 

(  \ 


Y  = 
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(5.9) 


and 


r  \ 
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dfK 

f  \ 

(L> 

Y  = 

I 

\  V 

=  Z —  = 

dY 

0 

!  / 

(5.10) 


where  the  three  frequencies,  a»„  correspond  to  the  rate  of  change  of  each  of  Vinti’s  action- 
angles  coordinates.  As  is  consistent  with  the  action-angle  formulation,  this  says  the  angles, 
6j,  increment  linearly  in  time  according  to  these  basis  frequencies  while  the  conjugate 
momenta,  2),  remain  constant. 

With  this  in  mind,  return  to 


6X  =  Z 


d2<H 

flp  «  +  Xpert 


(5.4) 
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and  change  variables 


d-<K  .  dY  . 

b Y t0ial  -  Z-^^SY  +  Ydrag  -  AdY  +  ^^Xjrag 


(5.11) 


where 


A  = 


0 


\ 

d(o 

dl 


0  0 


(5.12) 


v  ) 

The  Jacobian  matrix  relating  small  changes  in  the  Cartesian  frame  to  those  in  action-angle 
variables, 


SX  d  (r,  v) 
3Y  =  d(0, 1) 


(5.13) 


and  the  frequency  partials  matrix,  |y,  are  calculated  numerically  within  Wiesel’s  solution 
using  a  similar  method  to  the  one  presented  in  §3. 1.5.1  for  finding  the  Phi  matrix.  Then, 
assuming  this  mapping  between  X  and  Y  is  one-to-one,  continuous,  and  differentiable 
between  0  and  2n.  the  matrix  ||  is  inverted  to  change  drag  accelerations  into  action-angle 
form.  The  forcing  term  in  Equation  5.11  then  becomes 
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(5.14) 


It  follows  that 


dd\ 

£0drag  =  (  ^  |  Ud 


drag 


(5.15) 


and 


drag  —  ^  j  &drag  (5.16) 

We  now  have  an  expression  for  the  change  to  the  Vinti  action-angle  solution  due  to  air  drag. 
It  is  with  this  system  of  equations  that  an  analytical  solution  will  be  sought. 


5.4  Express  as  Fourier  Series 

The  expected  behavior  of  this  forcing  term  will  be  assumed  to  be  periodic  or,  at  the 
least,  contain  periodic  components.  This  is  reasonable  considering  the  local  changes  in 
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atmospheric  density  experienced  by  a  satellite  around  a  non-circular  orbit.  At  perigee, 
the  density  is  at  a  maximum  and  reduces  to  a  minimum  at  apogee.  This  value  oscillates 
between  minimum  and  maximum  as  the  satellite  revolves  at  a  period  commensurate  with 
that  of  the  orbital  period. 

Since  there  is  periodic  behavior  in  this  set  of  nonlinear,  coupled  equations,  it  should 
be  easily  expressed  as  multiple  Fourier  series  that  are  functions  of  6.  To  demonstrate 
finding  the  coefficients  needed  for  this  formulation,  consider  the  changes  per  B*  due  to 
drag  (denoted  by  a  prime), 

s  =  f(«)  (5.17) 


Changes  to  the  system  are  considered  on  a  per  B*  basis  for  number  conditioning 
reasons.  Return  to 


A] 


tag 


and  considering 
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1  CrfApref  p 

—  m  /  Vcf 

_  CdApref 
2m 


vv 


(5.6) 


(5.7) 


notice  that  the  reference  air  density  is  combined  into  the  leading  B*  parameter  and  then 
divided  out.  Considering  each  term  within 


adrag  =  -B*  —  |v|v  (5.8) 

Pref 

the  drag  coefficient  is  very  small  (B*  «  1),  yet  constant.  However,  atmospheric  density 
at  the  satellite’s  altitude  divided  by  the  reference  density  is  near  order  unity.  Likewise, 
the  term  |v|v  is  of  similar  order  because  dimensionless,  earth-based  distance  and  time 
units  result  in  a  magnitude  of  velocity  near  1  for  most  low  earth  orbits.  Thus,  it  is 
numerically  desirable  to  analyze  and  solve  the  function  Ydrag  divided  by  B*.  After  a  solution 
is  calculated  for  a  certain  size,  shape,  and  orientation  of  orbit,  this  constant  can  then  be 
multiplied  back  in  for  a  specific  satellite. 
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Looking  at  Equation  5.6  for  the  acceleration  due  to  drag,  it  does  not  appear  to  contain 
a  direct  relation  with  the  angle,  dt.  This  functional  dependence  will  be  computed  by 
incrementing  the  three  angles  from  0  to  2 n  in  succession,  calculating  velocity  and  density 
values  along  the  way  resulting  in  acceleration  quantities  for  each  angle  combination.  Next, 
the  following  equations 


Cj  = 


f(0)  *  cos(j  •  6)ddidd2d63 
f (6)  *  sin(j  •  6)ddidd2d62 


(5.18) 


are  solved  using  a  rapid  and  straightforward  Simpson’s  rule  integration.  An  elaborate 
integration  scheme  is  not  necessary  here  due  to  the  expected  slow-changing  behavior  of 
the  forcing  term  as  a  function  of  theta. 

After  the  coefficients  are  in  hand  and  following  Wiesel’s  notation[35:191]  (see  next 
section),  this  series  looks  like 

(  \  (  .  \ 
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01* drag 

d\  B" 


Co  +  ^  [CjCO.s'  (j  •  cot)  +  Sj  sin  (j  •  m/)] 

j*o 


(5.19) 


This  representation  is  a  set  of  six  Fourier  series  where  Cj,  Sj  are  vector  coefficients  with 
j  being  a  vector  summation  index1 1 .  The  dot  products  j  •  co  are  linear  combinations  of  the 
basis  frequencies.  These  frequencies  multiplied  by  time  give  the  angle  coordinates. 

In  this  form,  finding  a  solution  to  these  equations  is  straightforward.  Analytically 
integrate  Equation  5.19  with  respect  to  time  to  result  in  the  system  of  equations 
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(5.20) 


ii 


Due  to  sine/cosine  ambiguities  (sin  —x  =  -  sinx,  etc.),  the  first  non-zero  index  shall  not  be  negative 
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Moving  forward,  it  will  be  necessary  to  make  some  abbreviations  for  brevity.  Return 


to  Equation  5.19  and  these  key  elements  of  the  Fourier  series  may  be  expressed  as 
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with  the  corresponding  solution  being 
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(5.22) 

The  lead,  constant  coefficient  will  use  the  same  notation  because  it  does  not  change  when 
integrated. 

This  section  showed  how  the  forcing,  drag  accelerations  can  be  expressed  and 
manipulated.  Including  these  changes  in  the  global  solution  requires  an  extra,  deliberate 
step,  to  be  presented  in  §5.6. 


5.5  An  Example  Using  Vector  Indices 

An  example  of  using  summation  indices  in  this  fashion  will  now  be  offered  to  provide 
clarity.  Consider  one  of  the  six  equations  to  be  expressed  with  this  method.  For  instance, 
changes  to  6\  due  to  drag  could  resemble  a  multiply  periodic  relation  such  as 

0X  =2.4-4.1sin(-02  +  03)  +  O.3cos(01  -202)  +  1.7cos(0!  -  03)  (5.23) 


Table  5.1  illustrates  how  these  would  be  stored.  Since  0,-  =  ojjt,  the  integrated  solution  with 
respect  to  time  becomes 
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(jl)  1  —  2a>2 


sin  (0j  -  202)  + 


1.7 


a>\  -  a»3 


sin  (0!  -  03 ) 


(5.24) 


Table  5.2  shows  how  these  solution  indices  would  be  stored.  These  indices  can  now  be 
summed  whenever  a  solution  at  a  specific  time  is  desired  using  the  summation  convention 
presented  earlier. 
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Table  5.1:  Expressing  Summation  Quantities 
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Table  5.2:  Expressing  Solution  Summation  Quantities 
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5.6  Small  Changes  to  a  Coupled  System 

Now  that  the  forcing  term  is  expressed  by  a  multiple  Fourier  series,  changes  to  the 


system  as  a  whole  can  be  considered.  Return  to 
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Notice  that  the  angle  equations  are  coupled  with  those  of  momenta.  Luckily,  the  latter  can 

be  solved  directly  because  the  only  change  in  momenta  is  due  to  drag. 

Return  to  the  Fourier  solution 
/  \ 
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and  only  consider  the  momenta  series.  Using  abbreviations  introduced  in  §5.4,  this  solution 
can  be  expressed  as 


^t'otai  =  srdme  =  cut  -  to)  +  rsj  „  =  Clt  +  tsm  -  6i'(t0). 


(5.27) 


where  6I'drag(to)  combines  constant  and  periodic  contributions  to  momenta  at  the  initial 
time.  Then,  inserting  this  into  the  total  angles  solution.  Equation  5.26  becomes 

dco 
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(5.28) 

Separate  constant  from  periodic  terms  and  integrate  this  with  respect  to  time,  leaving 


<)ta,  =  ^  (^(r  -  t20)  -  6I\to)(t  -  f0))+Cg(f-f0)+ J  SAt)  +  TS^dt  (5.29) 

This  shows  that  a  linear  drift  in  the  momenta  from  the  lead,  constant  coefficient  becomes 
quadratic  in  the  angle. 

Now,  combine  like  coefficients  while  expanding  the  periodic  portions 
^t'otal  =  ^  { ^Cl(t2  -  t2)  +  6I'(to)(t  -  to))  +  Cg(r  -  to) 
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then  integrate  the  periodic  portion 
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Equations  5.27  and  5.31  complete  incorporating  small  changes  due  to  air  drag  into 
perturbations  of  the  global  solution.  These  small  changes  can  then  be  converted  back  to 
Cartesian  space  and  added  to  the  predicted  state  as  provided  by  the  Vinti  solution  at  a  given 
time.  This  conversion  is  performed  using  previously  supplied  matrix,  from  Equation 
5.13.  Incorporating  these  changes  into  orbit  fitting  as  explained  in  Chapter  III  will  be 
discussed  in  §5.9. 

5.7  Possible  Resonances 

In  practice,  for  all  orbits  attempted,  the  three  basis  frequencies  returned  by  the  Vinti 
solution  were  very  nearly  equal  each  other.  The  difference  between  any  two  are  on  the  order 
of  10-3.  Considering  the  summation  indices  that  account  for  all  the  linear  combinations  of 
frequencies,  there  are  certain  combinations  that  could  be  problematic. 

For  instance,  the  frequency  combination  that  corresponds  to  the  index  (1, 0,  -1)  would 
mean  <z>3  would  be  subtracted  from  to  i  resulting  in  a  very  small  number.  Looking  at 
the  momentum  solution  in  Equation  5.20,  this  number  is  in  the  denominator  of  the  new 
coefficients.  In  the  original  frequency  analysis,  if  any  of  the  coefficients  are  non-negligible 
this  new  term  will  become  noticeable.  This  effect  is  then  exacerbated  when  the  momenta 
solution  is  then  integrated  again  to  provide  the  angle  solution.  See  how  this  term  is  then 
squared  in  Equation  5.31.  Now,  a  term  that  was  merely  noticeable  before  is  now  significant. 
This  indicates  that  there  is  a  possibility  for  resonances  when  accounting  for  drag  using  this 
method. 

As  an  example,  consider  Orbit  2.  Table  5.3  shows  a  couple  sets  of  coefficients  found 
for  the  first  of  the  six  equations.  Anything  smaller  than  10-8  is  approximated  as  zero.  It 
is  important  to  note  that  the  majority  of  the  coefficient  pairs  are  zero  for  the  various  angle 
combinations. 
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Table  5.3:  Example  Coefficients  For  Orbit  2,  Coordinate  1 


j\ 

J2 

h 

Q 

3 

1 

2 

0 

0 

1.64e-2 

2 

-2 

0 

0 

2.01e-2 

Considering  the  resulting  basis  frequencies  for  this  orbit  are  (oT  =  (0.85042, 0.85 126, 0.85047), 
Table  5.4  shows  the  new  coefficients  after  analytical  integration  of  the  equations.  Notice 
how  for  the  first  angle  combination,  the  absolute  value  of  the  coefficient  decreased  by  more 
than  half  as  a  result  of  dividing  by  io\  +  2cl>2  or  (0.85042  +  2  *  0.85126).  However,  for  the 
second  combination,  it  is  divided  by  2oj\-2oj2  or  2  *(0.85042 -0.85 126)  =  -0.00168.  This 
increases  the  value  of  the  coefficient  by  four  orders  of  magnitude  from  2.01e-2  to  -5.86e+2. 


Table  5.4:  Example  Integrated  Coefficients  For  Orbit  2,  Coordinate  1 


j  i 

h 

h 

Cj 

1 

2 

0 

-6.90e-3 

0 

2 

-2 

0 

-5.86e+2 

0 

5.8  Atmosphere  Model 

This  investigation  uses  an  atmosphere  model  developed  by  Regan  and  Anandakrish- 
nan.  Instead  of  a  purely  exponential  model  for  density  as  a  function  of  reference  density 
and  height,  it  breaks  the  atmosphere  into  different  layers  of  strata.  Using  a  hybrid  method¬ 
ology,  it  uses  the  US  1976  Standard  Atmosphere  for  altitudes  from  0  to  86  km  and  the  US 
1962  Standard  Atmosphere  model  for  altitudes  above  86  km.  This  older  model  for  higher 
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altitudes  simplifies  calculations  greatly  by  using  linear  lapse  rates  instead  of  exponential 
and  elliptical  ones  found  in  two  of  the  higher  strata  levels.  These  minor  differences  have 
been  argued  by  Regan  and  Anandakrishnan  to  pale  in  comparison  to  other  changes  in  at¬ 
mospheric  density  that  are  less  predictable,  e.g.  incoming  solar  radiation.  Being  consistent 
with  the  motivation  of  the  current  research,  this  approach  was  chosen  to  maintain  simplis¬ 
tic,  rapid  computer  processing  while  demonstrating  a  relatively  accurate  solution  method. 
[78] 

In  reality,  the  altitude  used  to  determine  this  density  profile  pertains  to  height  above 
the  geoid.  Constant  fines  of  density  actually  take  on  the  shape  of  earth,  which  is  not 
spherical.  Therefore,  extra  calculations  are  needed  to  determine  local  altitude  relative  to  the 
ground  at  any  given  time.  The  current  research  will,  however,  only  use  spherical  altitude. 
Obviously,  there  may  exist  slight  differences  in  exact  densities  by  neglecting  nonspherical 
effects.  However,  the  truth  model  and  perturbation  routine  will  both  make  this  assumption. 
This  commonality  can  still  provide  a  validation  of  the  method  presented. 

5.9  Incorporating  Drag  Into  Orbit  Fitting 

Recall  from  Chapter  III  how  the  Phi  matrix  is  calculated  at  each  observation  then  used 
for  orbit  fitting.  Changes  to  the  six  state  variables  at  some  epoch  time  are  related  to  those 
at  an  observation  time  through  the  6x6  state  transition  matrix.  This  is  then  used  in  the 
accumulated  equation 

SX(t0)  =  (TrQ~l  T)-1  TtQT 1  r  (5.32) 

which  solves  for  six  corrections  to  the  state  variables  using  only  three  residuals  (remember 
T  =  H ®  from  Equation  3.7).  This  mechanism  allows  for  a  convenient  avenue  by  which  to 
calculate  and  store  the  necessary  quantities  needed  to  adjust  the  Vinti  orbit  as  affected  by 
drag. 

Add  a  seventh  state  to  solve  for  B* .  Now,  the  modified  state  vector  is 
X'  =  (x,y,z,x,y,z,B*)T  and  B  =  0.  Perform  partials  of  the  state  at  time  t  with  respect  to 
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the  state  at  to  to  get  the  modified  Phi  matrix 


dx  (t) 
dx(to) 

dx(t) 

dytto) 

dx(t ) 
dz(t0) 

dx(t) 
OB * 

dy(t) 

dxtto) 

dy(t) 

dytto) 

dytt) 

ditto) 

dytt) 
dB * 

ditt) 

dx(to) 

ditt ) 
dytto) 

ditt) 

ditto) 

ditt) 
OB * 

0 

0  ., 

..  0 

1 

(5.33) 


The  bottom  row  is  0  because  the  drag  coefficient  is  considered  a  constant.  Notice  the  new 
values  in  the  last  column.  These  are  changes  in  the  X  state  with  respect  to  B*.  However, 
the  formulations  in  previous  sections  solved  for  changes  to  the  action-angle  variables  per 
B*.  Arriving  at  these  A- space  quantities  requires  using  the  previously  found  partials  matrix 
(see  Equation  5.13)  to  perform  a  change  of  variables  from  theta  and  momenta  to  position 
and  velocity. 

The  additional  fitting  process  to  incorporate  drag  for  a  given  observation  follows. 
Before  the  fitting  routine  kicks  off,  perform  the  Fourier  analysis  discussed  earlier  to 
calculate  coefficients  for  the  six  series  as  in  Equation  5.18.  These  involve  constant  and 
periodic  coefficients.  Return  to  Equation  5.20  and  solve  each  of  the  six  equations  at  a  given 
time  with  the  drag  coefficient  divided  out.  This  solution  is  rapid  considering  only  time 
needs  to  be  supplied  for  multiplication  then  summing  over  the  Fourier  series.  These  values 
provide  incremental  quantities  of  the  action-angle  variables  per  B*.  Convert  these  values  to 
Cartesian  space  by 

d(r,\)(t)  _  <9(r,v)  6(6, 1)(t) 

dB*  ~  8(0,1)  B*  [  ’ 

and  save  them  in  the  Phi  matrix. 

In  parallel,  the  original  Phi  parameters  are  determined  normally.  The  fitting  scheme 
then  continues  for  each  iteration  as  described  earlier  until  convergence  is  declared.  B* 
is  now  a  “solve  for”  parameter  giving  the  routine  another  setting  to  adjust  in  hopes  of  an 
optimal  orbit  solution.  This  is  reasonable  considering  the  difficulty  involved  in  determining 


68 


a  craft’s  B*  value  at  a  specific  time.  This  value  is  in  theory  a  constant  but  realistically  as 
the  satellite  changes  orientation,  the  frontal  area  relative  to  its  direction  of  travel  can  vary. 
Further,  accurately  determining  a  vehicle’s  drag  coefficient  (Cp)  poses  difficult  without  a 
wind  tunnel  that  can  replicate  orbital  velocities  and  densities. 

5.10  Results 

Test  orbit  1  (parameters  can  be  found  in  Table  3.1)  is  analyzed  using  this  procedure. 
This  orbit  is  low  in  the  atmosphere  and  should  provide  a  baseline  with  which  to  demonstrate 
capability.  This  approach  calculates  changes  in  theta  and  momenta  per  B*  so  results  should 
be  able  to  scaled  for  various  realistic  drag  parameters. 

Recall  that  the  development  of  the  variational  action-angle  relations  in  §5.6  revealed 
the  angle  solution  is  dependent  upon,  or  coupled  with,  that  of  the  momenta.  So  looking 
at  momenta  perturbations  first,  Figure  5.1  shows  the  total  change  due  to  drag  over  one 
day.  Notice  the  varied  effects  across  the  three  momenta.  Momenta  1  is  near  constant  while 
momenta  2  and  3  change  linearly.  The  change  in  momentum  3  is  more  pronounced  and 
also  exhibits  a  periodic  oscillation  of  seemingly  constant  amplitude.  Through  the  coupling 
matrix  Jy),  this  noticeable  change  in  momenta  affects  all  three  angles.  As  predicted  in 
Equation  5.31,  Figure  5.2  shows  that  this  linear  growth  manifests  as  a  quadratic  6  drift. 
The  first  angle  experiences  a  periodic  affect  but  all  three  have  near  equal  secular  quadratic 
growth.  Note  that  the  behavior  of  angles  2  and  3  in  the  figure  are  identical. 
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Hours 

Figure  5.1:  Momenta  Perturbations 


To  approximate  model  error,  a  fourth-order  Hamming  numerical  integration  algorithm 
accounting  for  air  drag  will  serve  as  a  truth  model.  It  is  modified  slightly  to  call  upon  the 
Vinti  solution  package  in  order  to  output  action-angles  periodically.  The  geopotential  of 
this  model  has  also  been  adjusted  to  mimic  that  of  Vinti’s.  Therefore,  the  only  perturbations 
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compared  with  this  model  are  those  due  to  drag.  Since  the  integrated  trajectory  is  perturbed 
by  the  same  atmosphere  model,  these  instantaneous  values  represent  osculating  Vinti 
parameters  and  can  then  be  compared  to  the  analytically  predicted  values  at  commensurate 
time  steps. 

Figures  5.3  and  5.4  show  the  differences  in  momenta  and  angle  between  the  current 
approach  and  the  integrated  truth  model.  These  differences  are  presented  as  percentage 
errors  of  the  calculated  truth.  Momenta  2  and  3  exhibit  a  linear  increase  in  percentage  error 
and  equal  about  3.5%  after  one  day.  It  is  interesting  to  note  that  the  expected  quantity  of 
momentum  3  is  about  5  times  greater  than  momenta  2  but  the  percent  error  is  identical  in 
magnitude  and  behavior  with  the  exception  of  the  initial  transient.  It  is  important  to  note 
that  erratic  behavior  displayed  in  some  errors  near  t  =  0  are  a  result  of  a  small  divisor 
in  calculating  percent  error.  The  difference  in  perturbations  is  divided  by  the  magnitude  of 
integrated  pertubation  which  is  0  at  the  initial  time  and  very  small  for  the  first  few  timesteps. 
Figure  5.5  shows  error  growth  for  a  more  extended  time.  The  predicted  angle  perturbations 
are  approximately  87%  correct  after  5  days. 

5.11  Conclusions 

Consistent  with  the  motivations  of  the  current  research,  numerical  integration  of  the 
equations  of  motion  has  been  avoided.  What  resulted  is  a  straightforward  formulation 
of  how  air  drag  should  affect  the  reference  orbit.  This  solution  allows  for  a  relation  to 
be  determined  ahead  of  orbit  fitting  that  results  in  a  function  of  time.  Linear  changes 
in  momenta  resulted  in  quadratic  angle  as  expected.  The  angle  solution  from  techniques 
presented  in  this  chapter  remain  97.5%  and  87%  correct  after  1  and  5  days,  respectively. 
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Figure  5.4:  Error  in  Angle  Perturbations 
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VI.  Air  Drag  -  An  Alternate  Approach 


The  author’s  motivation  for  implementing  the  type  of  perturbation  technique  presented 
in  Chapter  V  was  having  the  action-angle  variables  available  to  generate  an  integrable 
analytical  function.  As  a  validation  of  this  technique  using  the  Vinti  solution,  an  alternate 
approach  is  presented  and  results  are  compared.  This  chapter  will  implement  a  similar 
formulation  using  another  dynamical  system  with  actions  and  angles. 


6.1  General  Perturbations  Using  Two  Body  Action- Angles 

As  discussed  earlier  in  §2.7,  obtaining  these  unique  action-angle  quantities  for  the  two 
body  problem  using  the  Hamilton-Jacobi  theory  results  in  Delauney  variables.  Although 
the  two  body  problem  is  fully  degenerate,  the  one  available  frequency  (due  to  the  mean 
anomaly)  is  the  most  relevant  to  air  drag  perturbations  since  the  dissipating  effect  is  mostly 
in-track  or  in  the  velocity  direction.  It  is  intuitive  to  imagine  the  satellite  would  slow  down 
in  its  orbital  plane  due  to  air  drag  more  than  the  node  would  be  shifted  for  an  inclined  orbit. 

Using  two  body  problem  variables  with  equations  of  motion 

(  \ 


X-TBP  ~  f(X,  t )  — 


V 


[  -  /jr/r3  j 

consider  air  drag  perturbations  as  a  forced  linear  system.  Linearize  around 
orbit  and  add  a  small  forcing  term  to  obtain  the  equations  for  small  changes  to 
Similar  to  before,  this  can  be  expressed  as 


(6.1) 

a  reference 
the  system. 


dX  —  A(t)6X  +  Xdrag 


where  A(t)  =  Vf  and 


A(t) 


0 

I 

^2,1 

0 

(6.2) 


(3.21) 
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For  the  two  body  problem,  the  gravity  acceleration  terms  found  in  the  lower  left  sub-matrix 
are 


A  2,1  - 


-4  + 


3/ja~ 


3/Jxy 

3  jixz 
r 5 


3 yUAV 

M  ,  W 

.A  I  S 


3 jjyz 

r5 


3 iixz 
E 

3  nyz 


.E  _i_ 

r3  r5 


(6.3) 


Another  way  to  express  the  variational  changes  of  the  physical  variables  is  by  relating 
them  to  Delauney  elements  through  a  matrix  of  partial  derivatives.  This  conveys  as 


6X  = 


dX 

8D 


dD 


(6.4) 


Remember  from  before  that  these  elements  are  D  =  (L,  G,  H,  /,  g,  h)  where  L,  G,  H  are 
constant  momenta  and  /,  g,  h  are  corresponding  angles  and  are  related  to  the  orbital  elements 
as 


L  =  y/JIa 

l  =  M 

G  =  LV 1  -e2 

g  =  O) 

(2.1) 

H  =  G  cos  i 

h  =  Q 

For  the  two  body  problem,  all  values  are  constant  except  the  mean  anomaly,  M,  which 
varies  according  to  mean  motion,  n  =  Now,  take  a  time  derivative  of  Equation  6.4, 


6X  = 


(  d  dX\  dX  . 


(6.5) 


Possessing  two  ways  to  express  the  same  quantity,  set  Equation  6.2  equal  to  Equation  6.5, 
substitute  for  SX, 


d  dX  _  dX  .  dX  _  . 

7tdD6D  +  dDSD  ~  M0dD6D  +  Xdmg 


(6.6) 


and  solve  for  the  element  rate  of  change, 

again  assuming  that  the  mapping  between  Cartesian  and  Delauney  variables  is  one-to-one, 
continuous,  and  differentiable  from  0  to  In.  Similar  to  before  with  Vinti,  we  now  have  an 
expression  for  the  change  to  the  TBP  action-angle  solution  due  to  air  drag  with  which  we 
can  seek  an  analytical  solution. 
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6.2  Another  Coupled  System 

At  first  glance,  this  expression  for  the  time  rate  of  change  of  small  changes  to  Delauney 
elements  as  perturbed  by  air  drag  does  not  seem  to  have  the  same  coupling  effect  on  the 
lead  term  as  the  expression  found  for  similar  changes  to  Vinti  action-angles.  From  Chapter 
V,  the  expression  for  total  variations  for  Vinti 

d2(K  .  dY  . 

dY  total  =  Z^SY  +  Ydrag  =  ASY  +  —  Xdrag  (5.11) 

shows  the  solution  is  multiplied  by  A  which  is 

/  \ 

0  % 

A  =  91  (5.12) 

0  0 

V 

and  constant.  This  alternate  approach  arriving  at  the  expression  for  TBP  elements  contains 
the  matrices  and  ATBP  which  are  actually  periodic  in  this  application.  However,  when 
subtracted  to  form  the  entire  leading  term  that  multiplies  dD  it  becomes  a  constant  matrix 
like  before  in  the  Vinti  relation.  However,  only  one  element  is  non-negligible  as  can  be 
seen  in  a  representative  output  of  this  matrix: 

2.27e-013  5.68e-014  8.53e-014  -2.65e+000  -1.63e-009  9.95e-014 

-1.99e-013  -8.53e-014  -2.84e-014  -1.75e-009  1.86e-009  7.11e-015 

-1.39e-017  5.55e-017  -1.10e-016  2.10e-013  -1.14e-013  -4.44e-016 

1.33e-015  l.lle-015  8.88e-016  -2.84e-014  1.14e-013  -2.78e-017 

1.33e-015  8.88e-016  8.88e-016  2.84e-014  1.42e-013  -2.78e-017 

8.88e-016  4.44e-016  8.88e-016  -5.68e-014  1.14e-013  l.lle-016 

This  turns  out  to  be  element  [1,4]  which  couples  the  solution  of  the  first  action  momenta 
into  changes  in  the  first  angle,  the  mean  anomaly.  This  is  unlike  before  in  Equation  5.12 
where  the  upper  right  block  of  the  matrix  was  non-zero  which  coupled  all  the  momenta 
solutions  back  into  all  angle  variations.  So,  all  small  variations  to  action  and  angles  for 
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Delauney  are  equal  only  to  the  effect  due  to  drag  with  the  exception  of  the  mean  anomaly 
and  this  is  coupled  with  the  corresponding  momenta.  This  becomes 


66\  —  cl  \  +  86\jrag 
\  — 


(6.8) 


and 


SI, 


total 


SI, 


drag 


(6.9) 


Return  to  Equation  6.7  and  compare  the  forcing  term  to  the  one  found  before  in 
Equation  5.11.  The  mechanics  of  the  air  drag  perturbation  are  the  same  but  differs  in 
the  relation  between  Cartesian  and  action-angle  coordinates.  We  are  now  ready  to  analyze 
these  expressions  and  calculate  a  solution  using  the  technique  presented  previously. 


6.3  Fourier  Series  Approach  Revisited 

Possessing  the  forcing  relation  between  Cartesian  and  Delauney  variables  due  to  air 
drag  as  well  as  the  coupling  effect  between  mean  anomaly  and  its  momenta,  analysis  will 
follow  the  steps  prescribed  in  §5.4.  Using  the  same  atmosphere  model  and  take  the  new 
forcing  term 


D 


drag 


B* 


=  D) 


drag 


m 


(6.10) 


and  express  it  as  a  Fourier  Series 


(  \ 
80' 


SI 


)  drag 


(  \ 
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81  Vdrag 

8\  B“ 


Co  +  ^  [Cjco.s'  (j  •  col)  +  Sjsin  (j  •  m/)] 

po 


(5.19) 


Solving  as  before  this  becomes 


t  \ 

SQ' 


SI' 


=  C0(t  -  h)  +  ^  |  H ~sin  0  ’ ~  rJr.cos  (J  ’  "0 


/drag 


UO  L 


j  m 


j  co 


(5.20) 


Since  five  of  the  terms  are  not  coupled,  this  is  their  direct  solution  for  changes  due  to 
air  drag.  In  the  case  of  mean  anomaly,  the  coupling  term  needs  to  be  considered.  With  the 
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momenta  solution  inserted  and  recalling  the  notation  presented  in  §5.4,  the  angle  time  rate 
of  change  is 

< total  =  Cl  (cJ0lt  +  rSTl(t)  -  6I\(t0))  +  cj  +  TSfh  (6.11) 

Integrating  this  analytically,  the  solution  becomes 


^ I  .total  =  ci 


-Cq‘t  +  -  to)j  +  Cj(f  -  to) 


+X 

j*o 


SI 

J  O) 


-  Cl 


STl 
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(j  -")2 


sin  (j  •  tot)  - 


S01 

_ j_ 

j  -to 


+  C 1 


Cf1 


(j  ’^)2 


COS  (j  •  Oil ) 


(6.12) 


after  expanding  periodic  terms  and  combining  like  terms. 

Equations  5.20  and  6.12  complete  the  solutions  for  small  changes  to  a  TBP  system 
under  the  influence  of  air  drag  only.  We  now  possess  a  method  to  calculate  perturbed 
Delauney  actions  and  angles  to  compare  with  those  from  Vinti. 


6.4  TBP  Results 

For  the  sake  of  comparison,  test  orbit  1  will  be  analyzed  again  using  the  current 
method.  Identical  orbits  examined  through  two  different  lenses  provide  a  basis  to  draw 
some  conclusions  on  performance  and  behavior.  Examining  momenta  first,  Figure  6.1 
shows  momenta  perturbations  after  one  day. 
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Momenta  1  - 

Momenta  2 
Momenta  3 


Figure  6.1:  TBP  Momenta  Perturbations 


Comparing  Figure  5.1  to  6.1,  the  noticeable  difference  from  what  was  presented  for 
Vinti  is  how  all  momenta  are  affected  now.  Earlier,  only  the  third  momenta  incremented 
linear  in  time  and  the  rest  were  nearly  constant.  Remembering  the  expressions  in  Equation 
2. 1  for  the  TBP  momenta  this  makes  sense.  Each  quantity  is  directly  proportionate  to  the 
semi-major  axis,  a.  The  first  momenta,  L,  is  yfjla.  The  second  momenta  is  multiplied 
by  L,  and  the  third  momenta  is  multiplied  by  the  second.  Therefore,  as  the  dissipative 
effect  of  air  drag  reduces  the  orbits  energy  and  ultimately  altitude,  each  momenta  decreases 
accordingly.  The  rate  of  decrease  will  differ  as  functions  of  eccentricity  or  inclination. 
Otherwise,  the  magnitude  of  the  effect  is  very  similar.  Differences  from  truth  will  be 
presented  shortly  then  compared. 

Now  consider  changes  to  the  angle.  Before  with  Vinti,  each  angle  experienced 
quadratic  effects  due  to  the  coupling  effect  with  the  one  linear  momenta  perturbation. 
Figure  6.2  shows  the  results  for  TBP  angles.  As  expected,  the  dominant  effect  manifests 
itself  in  the  mean  anomaly  or  parallel  to  the  direction  of  travel.  Again,  the  angle  exhibits  a 
quadratic  shift  with  a  periodic  behavior  consistent  with  the  orbital  frequency.  The  node  and 
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argument  of  perigee  are  mainly  constant  with  a  slight  periodic  oscillation  with  the  orbital 
frequency. 


Using  the  same  truth  model  as  before  but  altered  to  output  TBP  quantities,  results  are 
compared.  Figures  6.3  and  6.4  show  the  error  in  the  current  method  as  a  percentage  of  the 
total  calculated  truth  perturbation.  The  error  growth  compared  to  the  total  value  of  expected 
perturbations  seems  stable.  The  percentage  growth  of  the  momenta  is  linear  and  is  on  the 
order  of  3%  after  one  day. 
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Figure  6.3:  TBP  Momenta  Error 
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Figure  6.4:  Mean  Anomaly  Error 


As  with  the  momenta  (Figure  6.3),  the  mean  anomaly  percent  error  growth  (Figure 
6.4)  appears  linear.  However,  a  more  noticeable  periodic  behavior  is  present  again 
corresponding  with  the  orbital  frequency.  Figure  6.5  extends  this  prediction  out  five  days 
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for  perspective.  The  data  suggest  that  this  perturbation  method  remains  about  98%  accurate 
after  one  day  and  88%  accurate  after  five  days  with  seemingly  stable  behavior.  With  this 
method,  comparing  error  growth  in  the  second  and  third  angles  is  not  necessary  since  there 
is  no  secular  growth.  Further,  calculating  a  percent  error  is  not  feasible  since  it  would 
require  dividing  by  zero  when  the  truth  perturbation  vanishes  periodically. 
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Figure  6.5:  Mean  Anomaly  Error  Over  5  Days 


6.5  Conclusion 

Although  behavior  is  expectedly  different  in  general,  the  magnitude  of  total 
perturbations  expected  and  resulting  errors  are  strikingly  similar.  This  provides  evidence 
that  this  analytical  approach  consisting  of  Fourier  analysis  and  rapid  solution  calculation 
is  valid.  Furthermore,  the  twofold  validation  to  demonstrate  viability  proves  the  newly 
modified  Vinti  solution  is  well  suited  for  such  an  approach. 
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VII.  Conclusion 


7.1  Summary 

A  preliminary  examination  into  the  claimed  “other  solvable  solution”  has  been 
presented.  This  solution  is  viewed  as  a  candidate  for  a  more  comprehensive  orbit 
determination  tool  with  which  to  increase  accuracy  for  space  situational  awareness 
purposes.  In  Chapters  III  and  IV,  a  nonlinear  least  squares  routine  was  developed  then 
utilized  to  demonstrate  the  effectiveness  of  Vinti’s  modified  solution  for  orbit  fitting. 
Chapter  V  illustrated  a  novel  approach  for  performing  classical  perturbations  on  this  new 
solution  using  air  drag.  Chapter  VI  validated  this  approach  by  repeating  the  steps  on  a  well 
known  and  understood  system  with  similar  and  relevant  characteristics. 

Although  it  is  not  ready  to  be  inserted  into  a  workstation  at  JSPOC,  Vinti’s  solution  as 
modified  by  Wiesel  holds  great  promise.  Between  the  literature  review  into  SGP4  accuracy 
and  the  comparisons  made  within  this  research,  it  is  safe  to  say  that  orbit  fitting  using 
this  solution  returns  very  similar  results  to  what  is  gained  by  the  Air  Force’s  analytical 
propagator.  Of  course,  this  is  only  true  for  a  limited  orbital  regime. 

7.2  Recommendations 

Future  work  related  to  Wiesel’s  solution  has  a  bright  future  although  there  are  near- 
term  limitations  and  obstacles  to  overcome.  The  potential  capability  and  utility  greatly 
outweighs  these  initial  difficulties.  A  summary  of  efforts  required  to  tap  into  these  benefits 
follows. 

•  The  current  research  developed  and  utilized  a  numerically  calculated  state  transition 
matrix.  This  was  accomplished  due  to  the  ease  of  implementation  so  research  could 
begin  with  analyzing  Wiesel’s  solution  in  a  timely  fashion.  In  theory,  there  should  be 
minor  trade-offs  in  performance  due  to  computing  time  needed  to  simply  go  ahead 
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and  calculate  the  values  at  each  time  step.  In  practice  however,  there  is  a  significant 
different  between  using  SGP4  for  least  squares  fitting  (which  uses  an  analytically 
derived  Phi  matrix)  and  Vinti.  An  analytical  Phi  matrix  for  Vinti’s  solution  should  be 
explored  and  used  for  orbit  fitting. 

•  Performance  degradation  experienced  for  solutions  near  polar  orbit  should  be 
addressed.  This  inclination  range  is  popular  for  various  reasons.  In  order  for 
this  solution  to  be  widely  used  and  accepted  for  LEO,  the  inclination  limitation 
needs  removing.  It  is  the  author’s  belief  the  troubles  lie  in  Wiesel’s  numerical 
integrals  adversely  affected  by  the  negative  argument  of  perigee  rate  between  critical 
inclinations. 

•  The  current  research  has  only  used  the  solvable  solution  approximating  the  behavior 
induced  by  the  full  geopotential.  However,  a  perturbations  approach  to  incorporate 
more  of  these  effects  is  necessary.  Including  higher  order  geopotential  perturbations 
should  minimize  periodic  and  secular  growth  in  future  error.  This  would  result  in 
accurate  predictions  being  valid  for  an  extended  period  of  time. 

•  The  next  largest  contributor  to  secular  error  growth  over  time  for  low  earth  satellites 
is  air  drag.  Although  a  method  accounting  for  air  drag  and  a  notional  fitting  technique 
was  presented,  orbit  fitting  using  them  was  not  demonstrated.  This  should  be 
explored  and  could  possibly  be  undertaken  in  parallel  with  either  of  the  suggestions 
above. 

•  Accomplishing  the  three  previous  efforts  would  round  out  an  effective  and  accurate 
orbit  fitting  scheme  for  LEO  objects.  Research  should  then  move  to  incorporating 
other  perturbing  effects  necessary  for  higher  altitude  orbits.  These  include  luni-solar 
third  body  effects  initially.  Having  developed  a  method  to  account  for  the  sun  and 
moon  should  translate  easily  to  adding  other  desired  bodies  if  necessary.  Also,  the  air 
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drag  scheme  should  prove  adaptable  to  allowing  for  solar  wind/radiation  effects.  The 
Fourier  series  sampling  approach  would  only  need  to  sample  a  different  perturbing 
acceleration  expression. 

•  Another,  more  advanced  aspect  of  orbit  determination  could  be  examined  that  is 
unique  to  the  small  percentage  of  orbiting  objects  that  are  active,  maneuverable 
satellites.  The  ability  to  identify  and  model  a  maneuver  within  a  batch  of  observations 
would  significantly  increase  utility  for  using  Vinti  in  orbit  fitting.  Otherwise,  these 
batches  of  observations  would  lend  to  poor  results  until  an  separate  program  sifted 
through  the  observations  then  provided  Vinti  with  a  new  epoch  state  at  the  time  of 
maneuver.  That  would  only  increase  the  dependence  on  other  systems. 

•  When  all  the  efforts  are  complete,  a  well-rounded  and  capable  solution  should  be 
at  hand.  At  that  point,  a  more  realistic  and  operationally  relevant  process  should  be 
examined.  Computing  efficiencies  throughout  the  code  could  be  implemented  to  keep 
processing  time  to  a  minimum.  Also,  observation  relations  for  representative  SSN 
data  are  required.  Then,  actual  satellite  data  could  be  fit  for  true  system  performance. 

7.3  A  Parting  Thought 

With  the  exponential  increase  in  objects  orbiting  earth,  interest  in  inexpensive  orbit 
determination  coupled  with  accurate  conjunction  assessment  will  continue  to  rise.  This 
will  remain  a  top  priority  not  only  for  the  scientific  community  but  also  for  international 
leadership.  Considering  the  threat  to  US  satellites  by  orbital  debris,  USAF  General  Lynn  III 
said  “Without  space  systems,  many  of  our  most  important  military  advantages  evaporate” 
[79].  This  impetus  should  motivate  the  community  to  not  allow  analytical  theory  to  stagnate 
and  sit  backseat  to  that  of  numerical.  However,  the  marriage  of  analytical  (using  Vinti’s 
solution)  and  numerical  (using  Wiesel’s  modification  followed  by  Fourier  expansions) 
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methods  as  described  here  has  the  potential  to  meet  the  need  of  the  space  monitoring 
community. 


86 


Appendix  A:  Vinti’s  Hamiltonian 


Using  Vinti’s  coordinate  system,  the  conjugate  momenta  can  be  found  as 
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Appendix  B:  Extra  Figures 


Figure  B.l:  Vinti  Geopotential  Fit  Residuals:  Orbit  3 


Figure  B.2:  Vinti  Geopotential  Fit  Residuals:  Orbit  4 
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Figure  B.3:  Vinti  Geopotential  Fit  Residuals:  Orbit  5 


Figure  B.4:  Full  Geopotential  Fit  Residuals:  Orbit  3 
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Figure  B.5:  Full  Geopotential  Fit  Residuals:  Orbit  4 


Figure  B.6:  Full  Geopotential  Fit  Residuals:  Orbit  5 
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RMS  vs  Inclination  and  Eccentricity  at  1200  km  perigee  height 


Figure  B.7:  Orbit  Fit  Performance:  Perigee  Height  =  1200  km 


RMS  vs  Inclination  and  Eccentricity  at  2000  km  perigee  height 


Figure  B.8:  Orbit  Fit  Performance:  Perigee  Height  =  2000  km 
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Appendix  C:  The  Critical  Inclination 


Under  the  influence  of  J2  the  expression  for  the  secular  effect  on  the  argument  of 
perigee  is 
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and  setting  this  rate  equal  to  zero  becomes 
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which  means  this  occurs  at  inclinations 
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or  angles  of  63.4349°  and  116.5650°.  These  are  termed  the  critical  inclination  where  the 
argument  of  perigee  is  fixed  and  the  orbit  does  not  rotate  in  its  own  plane.  Below  and  above 
these  inclinations,  the  orbit  advances  or  the  rate  is  positive.  Between  the  two  inclinations 
and  straddling  polar  orbits,  the  argument  of  perigee  backs  up  or  the  rate  is  negative.  Figure 
C.l  qualitatively  illustrates  this. 


Figure  C.l:  J2  Effect  on  Apsidal  Rate 
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